VU5N/ 


MIT/WHOI 2010-08 


Massachusetts Institute of Technology 
Woods Hole Oceanographic Institution 







Joint Program 
in Oceanography/ 
Applied Ocean Science 
and Engineering 



1930 


DOCTORAL DISSERTATION 


Controls on Earthquake Rupture and Triggering 
Mechanisms in Subduction Zones 

by 

Andrea Lesley Llenos 


June 2010 





MIT/WHOI 

2010-08 


Controls on Earthquake Rupture and Triggering Mechanisms in Subduction Zones 


by 

Andrea Lesley Llenos 

Massachusetts Institute of Technology 
Cambridge, Massachusetts 02139 

and 

Woods Hole Oceanographic Institution 
Woods Hole, Massachusetts 02543 

June 2010 

DOCTORAL DISSERTATION 


Funding was provided by the National Science Foundation Division of Earth Sciences, National Defense 
Science and Engineering Graduate Fellowship, United States Geological Survey National Earthquake 
Hazards Reduction Program Award and the Woods Hole Oceanographic Hollister Research Fellowship 

and Academic Programs Office. 


Reproduction in whole or in part is permitted for any purpose of the United States Government. This 
thesis should be cited as: Andrea Lesley Llenos, 2010. Controls on Earthquake Rupture and Triggering 
Mechanisms in Subduction Zones. Ph.D. Thesis. MIT/WHOI, 2010-08. 


Approved for publication; distribution unlimited. 


Approved for Distribution: 

/lA^luMCt /4 > 

Maurice Tivey, 

Department of Geology and Geophysics 





Edward A. Boyle 
MIT Director of Joint Program 


James A. Yoder 

WHOI Dean of Graduate Studies 


CONTROI.S ON EARTHQUAKE RUFITJRE AND TRIGGERINCJ 
MECHANISMS IN SUBDUCTION ZONES 


By 

Andrea Ixisley Llenos 
Sc.B., Brown University. 2004 

Submitted in partial fulfillment of the requirements for the degree of 
Doctor of Philosophy 
at the 

MASSACHUSETTS INSTITUTE OE TECHNOLOGY 

and the 


WOODS HOLE OCEANOGRAPHIC INSTITUTION 
June 2010 


© 2010 Andrea L. Llenos 
All rights reserved. 

The author hereby grants to MIT and WHOl permission to reproduce and 
to distribute publicly paper and electronic copies of this thesis document 
in whole or in part in any medium now known or hereafter created. 


Author. 




Certified by...!fctd 


Joint Program in Oceanography/Applied Ocean Science and Engineering 
Massachusetts Institute of Technology and Woods Hole Oceanographic Institution 
11 J February 12,2010 

I r 


l.l. 



Jeffrey J. McGuire 

Associate Scientist, Department of Geology and Geophysics, WHOl 

Thesis supervisor 


Accepted by. 



Bradford H. Hager 


Professor, Department of Earth, Atmospheric, and Planetary Science, MIT 
Co-chair, Joint Committee for Geology and Geophysics 








2 



CONTROLS ON EARTHQUAKE RUPTURE AND TRIGGERING 
MECHANISMS IN SUBDUCTION ZONES 


By 

Andrea Lesley Llenos 

Submitted to the MITAVHOI Joint Program in Oceanography/Applied Ocean Science 
and Engineering on February 12, 2010 in partial fulfillment of the requirements for the 
degree of Doctor of Philosophy in Marine Geophysics 


Abstract 

Large earthquake rupture and triggering mechanisms that drive seismicity in 
subduction zones are investigated in this thesis using a combination of earthquake 
observations, statistical and physical modeling. A comparison of the rupture 
characteristics of M > 7.5 earthquakes with fore-arc geological structure suggests that 
long-lived frictional heterogeneities (asperities) are primary controls on the rupture extent 
of large earthquakes. To determine when and where stress is accumulating on the 
megathrust that could cause one of these asperities to rupture, this thesis develops a new 
method to invert earthquake catalogs to detect space-time variations in stressing rate. 
This algorithm is based on observations that strain transients due to aseismic processes 
such as fluid flow, slow slip, and afterslip trigger seismicity, often in the form of 
earthquake swarms. These swarms are modeled with two common approaches for 
investigating time-dependent driving mechanisms in earthquake catalogs: the stochastic 
Epidemic Type Aftershock Sequence model [Ogata, 1988] and the physically-based rate- 
state friction model [Dieterich, 1994]. These approaches are combined into a single 
model that accounts for both aftershock activity and variations in background seismicity 
rate due to aseismic processes, which is then implemented in a data assimilation 
algorithm to invert catalogs for space-time variations in stressing rate. The technique is 
evaluated with a synthetic test and applied to catalogs from the Salton Trough in southern 
California and the Hokkaido corner in northeastern Japan. The results demonstrate that 
the algorithm can successfully identify aseismic transients in a multi-decade earthquake 
catalog, and may also ultimately be useful for mapping spatial variations in frictional 
conditions on the plate interface. 
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Chapter 1: Introduction 

The largest earthquakes in the world occur on the megathrust of subduction zones, 
where almost 90% of the total seismic moment is released [Pacheco and Sykes, 1992], 
Therefore, understanding the controls on large earthquake rupture and the triggering 
mechanisms that affect earthquake occurrence, in general, is critical in order to accurately 
assess both long-term and short-term seismic hazards in these regions. A wide spectrum 
of deformation occurs in subduction zones, ranging from the seismic (e.g., low frequency 
events, seismic tremor, microearthquakes, moderate-to-great interplate earthquakes) to 
the aseismic (e.g., stable slip, slow slip events, afterslip). All of these processes alter the 
regional stress state, but it is often difficult to monitor stress accumulation in subduction 
zones, as much of the seismogenic zone is located offshore. 

The asperity model is commonly used to explain variations in large earthquake 
rupture and seismogenic behavior [Lay and Kanamori, 1981]. In this model, asperities, 
which may be related to frictional properties along the plate interface, lock during the 
interseismic period and then suddenly release the accumulated strain in an earthquake. 
This model predicts that long-lived (timescales of -millions of years) frictional 
heterogeneity along the thrust controls rupture behavior. Recent studies by Song and 
Simons [2003] and Wells et al. [2003] based on correlating global and historical 
earthquake data with fore-arc geologic structure illuminated in the gravity field have 
provided evidence to support this model, suggesting that along-strike variations in 
frictional properties are important first-order controls on the rupture of great subduction 
zone earthquakes. 

An alternative hypothesis is the seismic gap model [e.g., Thatcher, 1989], in 
which the occurrence of large earthquakes is controlled by short-lived (timescales of 
-lOOs of years) time-dependent stress heterogeneities. In this model, parts of the fault 
that have not ruptured in a long time have accumulated more stress and therefore are 
closer to failure than parts of the fault that have experienced recent earthquakes. Thus, 
large earthquakes tend to fill in the gaps where stress has not been released recently, and 
so rarely occur in the same place as the previous large earthquake. 
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This thesis investigates the controls on large earthquake rupture and explores the 
triggering mechanisms that drive seismicity in subduction zones. Important questions 
that this thesis addresses include: 

1) What factors influence where the largest earthquakes occur? 

2) W hat controls how large these ruptures can grow? 

3) How can we map spatial variations in frictional conditions on the plate interface 
that could potentially affect earthquake rupture? 

4) H ow can we detect stress changes that could trigger these earthquakes, particularly 
offshore where land-based geodetic resolution is low? 

To address the last two questions in particular, I develop a new technique that uses 
earthquake catalogs to identify space-time windows in which stressing rate changes due 
to transient aseismic processes are occurring. This novel way of producing estimates of 
stressing rate variations in space and time from seismicity data can be used in tectonic 
settings besides subduction zones and has other potential applications besides transient 
detection, including investigations of the physical processes that trigger earthquakes and 
improvements to short-term and real-time hazard assessment. 

Chapter 2 sets the overall framework for this investigation by providing evidence 
to support the asperity model of large earthquake rupture. Following up on the 
hypothesis that gravity lows can be a proxy for seismogenic behavior [Song and Simons, 
2003; Wells et al, 2003], I use a combination of M >7.5 earthquake observations and 
gravity anomaly data in subduction zones to examine the relationship between fore-arc 
geological structure revealed in the gravity field and the rupture extent of large 
subduction zone earthquakes. I demonstrate that large ruptures tend to stop in regions 
with positive gravity gradients by estimating a characteristic rupture length and 
directivity for each earthquake and comparing them with the local gravity field. This 
suggests that local increases in the gravity field can be related to physical conditions on 
the plate interface that favor rupture termination, such as a transition from velocity¬ 
weakening (frictionally unstable) to velocity-strengthening (frictionally stable) behavior. 
As the gravity anomalies reflect geologic structure such as fore-arc basins that have 
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formed over timescales on the order of millions of years [Fuller et al, 2006], this 
provides evidence that long-lived frictional heterogeneities (i.e., asperities) are 
responsible for controlling the rupture extent of large earthquakes. 

This result raises several questions pertinent to estimating the seismic hazard in a 
subduction zone, including how to identify and map these frictional conditions on the 
plate interface, and how to determine when and where stress is accumulating that could 
cause an asperity to rupture. Transient aseismic deformation, such as slow slip events, 
fluid flow, or afterslip, may also alter the regional stress state, and their occurrence 
suggests the presence of velocity-strengthening conditions [e.g., Miyazaki et al, 2004]. 
Geodetic observations provide a way to monitor changes in deformation and estimate the 
degree of seismic coupling [e.g., Nishlmura et al, 2004], but a large part of the 
seismogenic zone of a subduction zone is located offshore, where land-based geodetic 
instruments have little resolution and it is challenging to place seafloor instruments. 

An alternative is to instead monitor spatial and temporal changes in seismicity 
rate. Strain transients due to aseismic processes, such as fluid flow [e.g., Hainzl and 
Ogata, 2005; Bourouis and Bernard, 2007], slow slip events [e.g., Segall et al, 2006; 
Ozawa et al, 2007; Lohman and McGuire, 2007; Delahaye et al, 2009], and afterslip 
[e.g., Matsubara et al, 2005; Hsu et al, 2006] have been shown to trigger variations in 
regional earthquake rates. Thus, seismicity rate variations in earthquake catalogs can 
potentially be used as a proxy to detect stressing rate variations in regions with limited 
geodetic data. Chapters 3-4 of my thesis are devoted to developing, testing, and applying 
this idea in a method to invert seismicity catalogs for stressing rate variations in space 
and time. Finally, in Chapter 5, I demonstrate how this new method can be applied to a 
subduction zone to detect aseismic strain transients and identify spatial variations in 
frictional conditions. 

Chapter 3 begins by investigating the efficacy of two widely used approaches, the 
stochastic Epidemic Type Aftershock Sequence (ETAS) model [Ogata, 1988] and the 
physically-based rate- and state-dependent friction model [Dieterich, 1994], for modeling 
earthquake swarms, which are often triggered by aseismic transients. The ETAS model 
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uses empirical aftershock triggering laws to model earthquake occurrence as a point 
process that can be described with optimizable parameters such as background seismicity 
rate and aftershock productivity, but it lacks a quantitative way of relating seismicity rate 
variations to stressing rate variations. The rate-state model, on the other hand, has been 
used to map seismicity rate variations to stressing rate variations [Dieterich et al, 2000; 
Toda et al, 2002], but the stressing rate changes estimated from this model will reflect a 
combination of variations due to the underlying aseismic process, as well as those due to 
aftershock sequences. In this chapter, I show that the rate-state model predicts a 
relationship between aftershock productivity and stressing rate that is not observed in real 
data. I also demonstrate that earthquake swarms in various tectonic settings appear as 
anomalies relative to the ETAS model [Ogata, 1988], because the heightened stressing 
rate during the swarms causes a significant increase in background seismicity rate while 
the other aftershock parameters remain relatively constant. These observations enable us 
to specify a combined ETA S/rate-state model of seismicity rate that models both 
aftershock activity as well as variations in background seismicity rate due to aseismic 
processes and provides a direct relationship to stressing rate. 

In Chapter 4, I implement the seismicity rate model specified in Chapter 3 into a 
data assimilation algorithm to invert seismicity catalogs for estimates of space-time 
variations in stressing rate. I set up a state-space model that describes the system using an 
underlying state vector, consisting of background stressing rate, aseismic stressing rate, 
and rate-state variable ^that evolves over time. I then use an extended Kalman filter to 
estimate the time history of the state variables in a given number of spatial boxes. I 
evaluate the algorithm using a synthetic catalog generated with known stressing rate 
histories including an aseismic transient, and show that it can successfully detect when 
and where the transient occurs. I then apply it to an earthquake catalog from the Salton 
Trough in southern California, where a number of aseismic transients such as afterslip 
and shallow creep events have been geodetically detected [e.g., Williams andMagistrale, 
1989; Lohman and McGuire, 2007; Wei et al, 2009]. The algorithm successfully detects 
the largest geodetically-observed transient in the catalog (the 2005 Obsidian Buttes 
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transient), and the filter estimate of the peak stressing rate during the transient is within a 
factor of 5 of the estimate from a geodetically-derived slip model [Lohman and McGuire, 
2007], This demonstrates that this approach can successfully identify space-time 
windows in which aseismic transients occurred from a multi-decade earthquake catalog. 

Finally, in Chapter 5,1 return to subduction zones and apply this new seismicity- 
based transient detection method to a catalog from northeastern Japan to identify space- 
time windows where aseismic processes such as afterslip may be occurring. Afterslip 
was observed geodetically following the 4 major interplate thrust events that occurred in 
this catalog (1989 M 7.1, 1992 M6.9, 1994 M7.6, and 2003 M8.0) [e.g., Miura et ah, 
1993; Kawasaki et ah, 1995; Heki et al, 1997; Miyazaki et al, 2004]. I show that 
seismicity rate anomalies relative to ETAS following these events can be detected from 
the earthquake catalog alone. Several smaller seismicity rate anomalies are also detected 
that can be associated with postseismic slip following M6.3-6.5 earthquakes and 
precursory slip prior to the 1994 M 7.6 Sanriku-oki earthquake. These transients were not 
observed geodetically but can be corroborated with repeating earthquake analyses 
[Uchida et al, 2003, 2004]. Moreover, analysis of the 2003 Tokachi-oki earthquake 
indicates that this method may be able to distinguish between velocity-weakening and 
velocity-strengthening patches on the fault. Aseismic slip can be associated with 
frictionally stable, velocity-strengthening behavior [e.g., Miyazaki et al, 2004], and the 
filter correctly identifies the spatial box where the peak afterslip occurred as opposed to 
where the coseismic rupture occurred (indicating frictionally unstable, velocity¬ 
weakening behavior). These results indicate that with improvements in spatial resolution 
and offshore earthquake detection levels, this method can help map the frictional 
conditions on the plate interface that may control large earthquake ruptures, as well as 
enhance our ability to detect where and how stress is accumulating on the megathrust, 
especially in regions further offshore where geodetic resolution is limited. 
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Chapter 2: 1 nfluence of fore-arc structure on the extent of great subduction 
earthquakes 


Abdract 

Structural features associated with fore-arc basins appear to strongly influence the 
rupture processes of large subduction zone earthquakes. Recent studies demonstrated that 
a significant percentage of the global seismic moment release on subduction zone thrust 
faults is concentrated beneath the gravity lows resulting from fore-arc basins. To better 
determine the nature of this correlation and to examine its effect on rupture directivity 
and termination, we estimated the rupture areas of a set of M w 7.5-8.7 earthquakes that 
occurred in circum-Pacific subduction zones. We compare synthetic and observed 
seismograms by measuring frequency-dependent amplitude and arrival time differences 
of the first orbit Rayleigh waves. At low frequencies, the amplitude anomalies primarily 
result from the spatial and temporal extent of the rupture. We then invert the amplitude 
and arrival time measurements to estimate the second moments of the slip distribution 
which describe the rupture length, width, duration, and propagation velocity of each 
earthquake. Comparing the rupture areas to the trench-parallel gravity anomaly (TPGA) 
above each rupture, we find that in 11 of the 15 events considered in this study the TPGA 
increases between the centroid and the limits of the rupture. Thus local increases in 
TPGA appear to be related to the physical conditions along the plate interface that favor 
rupture termination. Owing to the inherently long timescales required for fore-arc basin 
formation, the correlation between the TPGA field and rupture termination regions 
indicates that long-lived material heterogeneity rather than short timescale stress 
heterogeneities are responsible for arresting most great subduction zone ruptures. 
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[i] Structural features associated with fore-arc basins appear to strongly influence the 
rupture processes of large subduction zone earthquakes. Recent studies demonstrated that 
a significant percentage of the global seismic moment release on subduction zone 
thrust faults is concentrated beneath the gravity lows resulting from fore-arc basins. To 
better determine the nature of this correlation and to examine its effect on rupture 
directivity and termination, we estimated the rupture areas of a set of 7.5-8.7 
earthquakes that occurred in circum-Pacific subduction zones. We compare synthetic and 
observed seismograms by measuring frequency-dependent amplitude and arrival time 
differences of the first orbit Rayleigh waves. At low frequencies, the amplitude anomalies 
primarily result from the spatial and temporal extent of the rupture. We then invert the 
amplitude and arrival time measurements to estimate the second moments of the slip 
distribution which describe the rupture length, width, duration, and propagation velocity 
of each earthquake. Comparing the mpture areas to the trench-parallel gravity anomaly 
(TPGA) above each mpture, we find that in 11 of the 15 events considered in this study 
the TPGA increases between the centroid and the limits of the mpture. Thus local 
increases in TPGA appear to be related to the physical conditions along the plate interface 
that favor mpture termination. Owing to the inherently long timescales required for fore¬ 
arc basin formation, the correlation between the TPGA field and mpture termination 
regions indicates that long-lived material heterogeneity rather than short timescale stress 
heterogeneities are responsible for arresting most great subduction zone mptures. 
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1. Introduction 

[ 2 ] The largest earthquakes in the world occur in subduc¬ 
tion zones, where almost 90% of the total global seismic 
moment is released [Lay and Bilek, 2007]. However, the 
amount of seismic coupling widely varies from one subduc¬ 
tion zone to another, ranging from predominantly aseismic 
subduction in the Marianas to seismic subduction in Alaska 
[e.g., Scholz and Campos, 1995; Lay and Bilek, 2007]. 
Understanding this variation in earthquake occurrence in 
circum-Paciftc subduction zones has been the subject of 
numerous studies. The correlations between earthquake 
occurrence and such factors as age of subducting oceanic 
lithosphere, amount of sediment, bathymetric features on 
the subducting slab, and convergence rate have been inves¬ 
tigated [e.g., Mogi, 1969; Kelleher and McCann, 1976; Ruff 
and Kanamori, 1980; Pacheco et al, 1993; Scholz and 
Campos, 1995; Abercrombie et al, 2001]. However, wide 
variability in seismogenic behavior exists not only between 
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different subduction zones but within individual subduction 
zones themselves. 

[3] The asperity model is commonly used to describe this 
variability in seismic behavior [Lay and Kanamori, 1981]. 
The moment release during great earthquakes is nonuni¬ 
form, and the areas of high moment release are known as 
asperities. These asperities may occur because of variability 
in frictional properties on the plate interface, which may 
lock during the interseismic period and suddenly release slip 
in an earthquake [e.g.. Lay et al., 1982]. An alternative view 
is that time-dependent stress heterogeneity is the dominant 
factor controlling the extent of great earthquakes. Numerical 
simulations demonstrated that even a fault with uniform 
frictional properties can generate a complex sequence of 
events that rupture different portions of the fault in each 
rupture rather than repeatedly at a fixed asperity [Shaw, 
2000]. In this model, large earthquakes preferentially 
nucleate at the edge of a previous large rupture and 
propagate in the opposite direction providing a natural 
explanation for the observation that large subduction zone 
ruptures are predominantly unilateral along strike [McGuire 
et al, 2002]. Moreover, Thatcher [1989] used historical 
estimates of rupture area for great subduction zone earth¬ 
quakes to argue that these events are rarely repeats of the 
previous big earthquake and instead fill in regions where 
stress accumulation has not been relieved recently (the 
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X (km) 



Final slip (m) 

Figure 1. Characteristic rupture ellipse for the 2000 
Tottori earthquake (solid black line), rupture directivity 
(arrow), centroid location (triangle), and hypocenter loca¬ 
tion (star) plotted on top of slip model of Iwata et al. [2000], 


seismic gap hypothesis). While both “quenched heteroge¬ 
neity” (the asperity model) and “dynamic heterogeneity” 
(the seismic gap hypothesis) likely influence the details of 
individual ruptures to some extent, it is very important to 
determine which is the dominant behavior in subduction 
zones as they have very different implications for the long¬ 
term seismic hazard at a particular location. 


[ 4 ] In a given subduction zone, the trench-normal con¬ 
trols on earthquake rupture are relatively well understood. 
The fault width is constrained by the updip and downdip 
limits of the seismogenic zone, at depths of 5-10 km and 
25-55 km, respectively [e.g., Byrne et al, 1988; Pacheco et 
al., 1993; Tichelaar and Ruff, 1993; Oleskevich etal., 1999; 
Lay and Bilek, 2007]. These limits result from transitions 
from velocity-strengthening to velocity-weakening behavior 
[Scholz, 2002]. These transitions likely occur because of 
changes in properties such as sediment strength and mineral 
composition due to changing pressure and temperature 
conditions [e.g., Byrne et al, 1988; Hyndman and Wang, 
1993; Oleskevich et al., 1999]. 

[5] However, subduction zone earthquakes can release 
large amounts of seismic moment because extremely long 
fault lengths are possible along a subduction zone. An 
important question then becomes: what controls the along- 
strike limits of these great earthquakes? What stops a 

7.0 earthquake from continuing to rupture along a 
subduction zone and becoming a 9.0 earthquake? Some 
studies have shown that transverse structures such as ridges 
or seamounts in the subducting lithosphere often fragment 
the subduction zone and may provide natural barriers to 
rupture [e.g., Mogi, 1969; Kelleher and McCann, 1976; 
Cloos, 1992; Kodaira et al, 2000, 2002]. However, such 
structures can also prove to be asperities, as Abercrombie et 
al. [2001] found in the case of the 1994 Java tsunamigenic 
earthquake, where the majority of the seismic moment 
release was centered on a subducted seamount. Therefore 
the relationship between subducting geological structure 
and the extent of large individual ruptures remains unclear. 


1 


■Q 0.5 


0 

0 90 180 270 360 0 90 180 270 360 

10-12 mHz 
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0.5 


0 

0 90 180 270 360 0 90 180 270 360 
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Figure 2. Amplitude measurements (square) and model (circle) from a synthetic line source test at 
different frequency bands. The 3-D point source synthetics were used to simulate a 70 km long unilateral 
rupture propagating to the east with a velocity of 2.5 km/s. Low-amplitude measurements at azimuths of 
270 confirm that rupture propagated away from the west. The amplitude anomaly increases with 
frequency. In general, the model fits the measurements well at lower frequencies, but at the 14-16 mHz 
band they begin to differ significantly. 
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Figure 3. Location and focal mechanisms of the 15 > 

7.5 events in data set. 


[6] The upper plate also plays an important role in 
earthquake rupture. McCaffrey [1993] emphasized the im¬ 
portance of fore-arc rheology in seismogenic behavior; 
strong fore arcs produce more large earthquakes because 
of their ability to store elastic strain energy. Recent studies 
by Song and Simons [2003] and Wells et al. [2003] have 
demonstrated that large subduction zone earthquakes occur 
preferentially in areas along the plate interface which are 
overlain by fore-arc basins. Song and Simons [2003] found 
that 80% of the cumulative seismic moment release in the 
20th century occurred in the 30% of the total area of 
subduction zone which exhibited strongly negative gravity 
anomalies, indicators of the presence of fore-arc basins. 
Fore-arc basins tend to form in strong, stable wedges and 
therefore reflect the mechanical and frictional properties 
along the plate interface over which they lay [Byrne et al ., 
1988; Fuller et al, 2006; Wang and Flu, 2006]. To what 
extent then can fore-arc structure influence the rupture of 
individual earthquakes? 

[7] This study aims to investigate the relationship 
between earthquake rupture propagation and fore-arc struc¬ 
ture in greater detail. Where do large events start and stop 
with respect to along-strike structures in the gravity field? 
To address tliis question, we estimated the finite source prop¬ 
erties of a set of 15 subduction zone thrust events > 7.5) 



and compared them to fore-arc structure revealed by maps 
of trench-parallel gravity anomalies (TPGA) constructed by 
Song and Simons [2003]. We find that as a rupture 
approaches its eventual extent, the TPGA increases. This 
correlation, which reinforces the observations of Song and 
Simons [2003], suggests that the stress and frictional heter¬ 
ogeneities along the plate interface that control the rupture 
of large subduction zone earthquakes are expressed in the 
fore-arc gravity field. 

2. Methodology 

[s] To evaluate the relationships between fore-arc struc¬ 
ture and individual ruptures, we require well constrained 
information on the spatial extent of rupture that can be 
determined in subduction zones worldwide. This is often a 
difficult observational problem because most subduction 
zone ruptures occur offshore and often with limited geodetic 
and near-field seismic data. Even for the best data sets, 
detailed finite fault inversions of great subduction zone 
earthquakes are highly sensitive to model parameterization 
and station coverage. An example of this is the 2003 8.3 

Tokachi-oki earthquake. Despite the abundance of quality 
seismic and geodetic data recorded during the event, by far 
the best data set ever for a M^, 8 subduction zone rupture, 
the finite fault models produced following the earthquake 
differ in characteristics such as number of asperities and 
orientation of rupture [Miyazaki et al, 2004a]. Some studies 
found that the rupture was oriented more along strike of the 
trench [Yamanaka and Kikuchi, 2003] while others found 
rupture areas that were half the size and oriented downdip 
[Yagi, 2004; Honda et al, 2004; Koketsu et al, 2004; 
Miyazaki et al, 2004a]. The differences between the slip 
models highlight the sensitivity of the results to the different 
fault parameterizations, constraints, and data sets used in 
each study. Even for the best combined seismic and geo¬ 
detic data sets, the rupture area is only constrained to within 
a factor of two owing to the limited offshore coverage. 
Moreover, body wave based studies often have poor con¬ 
straints on the seismic moment and slip distribution owing 
to their relatively high frequency band [Pritchard et al, 
2007]. 

[9] An alternative approach is to utilize seismic surface 
waves to constrain only the gross features of the rupture. 


Table 1. Characteristic Rupture Dimensions of the Events in This Study* 


Event 

Location 


L^., km 

'Tc, S 

VqI, km/s 

Directivity Ratio 

19920902 

Nicaragua 

7.6 

74 ± 7 

40 ± 1 

1.8 ± 0.1 

0.97 ± 0.10 

19940602 

Java, Indonesia 

7.8 

86 ± 7 

16 ± 2 

4.2 ± 1.3 

0.78 ± 0.19 

19941228 

Sanriku-oki, Honshu, Japan 

7.7 

161 ± 14 

20 ± 3 

6.1 ± 1.9 

0.77 ± 0.16 

19950730 

Antofagasta, Chile 

8.0 

121 ± 10 

33 ± 1 

1.3 ± 0.1 

0.35 ± 0.04 

19950914 

Copala, Mexico 

7.3 

74 ± 4 

14 ± 1 

3.2 ± 0.3 

0.62 ± 0.06 

19951009 

Jalisco, Mexico 

8.0 

77 ± 6 

32 ± 1 

2.1 ± 0.1 

0.86 ± 0.08 

19951203 

Kurile Islands, Russia 

7.9 

121 ± 9 

36 ± 1 

1.1 ± 0.1 

0.32 ± 0.03 

19971205 

Kamchatka, Russia 

7.8 

58 ± 2 

34 ± 1 

1.6 ± 0.1 

0.95 ± 0.02 

20010623 

Arequipa, Peru 

8.4 

167 ± 2 

27 ± 1 

1.1 ± 0.2 

0.18 ± 0.02 

20030122 

Colima, Mexico 

7.5 

92 ± 6 

29 ± 1 

1.4 ± 0.1 

0.44 ± 0.03 

20030925 

Tokachi-oki, Hokkaido, Japan 

8.3 

40 ± 2 

45 ± 1 

0.9 ± 0.1 

0.96 ± 0.02 

20031117 

Rat Islands, Alaska 

7.8 

70 ± 21 

29 ± 1 

2.2 ± 0.4 

0.91 ± 0.24 

20050328 

Sumatra, Indonesia 

8.7 

137 ± 19 

47 ± 3 

2.3 ± 0.3 

0.79 ±0.15 

20060717 

Java, Indonesia 

7.7 

108 ± 5 

55 ± 1 

1.9 ± 0.1 

0.98 ± 0.05 

20061115 

Kurile Islands, Russia 

8.3 

93 ± 3 

52 ± 1 

1.8 ± 0.1 

0.98 ± 0.02 


**Event number is date as year, month, day. Errors are ±1 standard deviation. 
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Table 2. Second Moments Inversion Results^ 


Earthquake 

^(0.2) 


AO-i) 





AP.O) 



r 

B 

<P 

rr 

r0 

np 

ee 

9(p 

(jxp 

Nicaragua 

406 ± 23 

-11 ± 3 

374 ± 19 

615 ± 14 

6 ± 1 

-10 ± 6 

-18 ± 7 

423 ± 225 

586 ± 131 

997 ± 135 

Java 

65 ± 19 

26 ± 4 

-70 ± 20 

262 ± 19 

56 ± 1 

214 ±46 

60 ± 17 

1514 ± 291 

-467 ± 132 

1140 ± 170 

Sanriku-oki 

103 ± 29 

4 ±2 

-292 ± 14 

-555 ± 18 

0.3 ± 0.3 

-4 ± 3 

-26 ± 9 

1093 ± 231 

990 ± 127 

6260 ± 1089 

Antofagasta 

280 ± 24 

22 ± 2 

351 ± 14 

-53 ± 17 

4.7 ± 0.4 

48 ± 5 

-26 ± 3 

1983 ± 272 

-1682 ± 271 

1963 ± 381 

Copala 

51 ± 5 

-10 ± 1 

-145 ± 4 

79 ±4 

3.1 ± 0.4 

18 ±4 

-44 ± 5 

601 ± 63 

237 ± 43 

1311 ± 126 

Jalisco 

259 ± 10 

12 ± 3 

-379 ± 14 

-373 ± 11 

8 ± 1 

-15 ± 10 

-18 ± 9 

842 ± 135 

686 ± 107 

730 ± 109 

Kurile 

320 ± 18 

-28 ± 4 

-306 ± 13 

164 ± 16 

4 ± 1 

29 ± 2 

-46 ± 5 

433 ± 124 

-434 ± 85 

3607 ± 563 

Kamchatka 

289 ± 6 

-41 ± 2 

310 ± 17 

-345 ± 16 

91 ± 13 

9 ± 8 

93 ± 8 

459 ± 42 

-363 ± 31 

486 ± 42 

Arequipa 

187 ± 16 

45 ± 5 

81 ± 11 

187 ± 22 

206 ± 2 

329 ± 16 

192 ± 45 

739 ± 14 

-737 ± 28 

6859 ± 177 

Colima 

210 ± 8 

0 ± 4 

-264 ± 12 

129 ± 8 

2.0 ± 0.4 

-7 ± 5 

-6 ± 3 

1726 ± 229 

717 ± 119 

725 ± 102 

Tokachi-oki 

500 ± 10 

-2 ± 3 

-434 ± 20 

-44 ± 8 

0.4 ± 0.2 

2 ± 2 

-0.4 ± 0.3 

406 ± 30 

26 ± 12 

178 ± 32 

Rat Islands 

208 ± 18 

-6 ± 3 

-59 ± 9 

-459 ± 55 

0.8 ± 0.1 

1 ± 3 

12 ± 3 

151 ± 208 

-15 ± 305 

1228 ± 725 

Sumatra 

559 ± 73 

15 ± 5 

1200 ± 42 

416 ± 35 

2 ± 1 

31 ± 11 

14 ± 5 

4503 ± 1031 

-660 ± 581 

2117 ± 854 

Java 

761 ± 31 

49 ± 5 

1022 ± 22 

1053 ± 24 

5 ± 1 

65 ± 5 

68 ± 7 

1435 ± 112 

1402 ± 61 

1613 ± 321 

Kurile 

667 ± 12 

-50 ± 5 

-1077 ± 38 

472 ± 20 

41 ± 1 

71 ± 7 

-31 ± 3 

1804 ± 106 

-788 ± 50 

431 ±26 


‘^Errors are ±1 standard deviation. Units of are s^, units of are km s, and units of are km^. 


such as its extent and directivity. Surface waves have the 
natural advantages of complete azimuthal coverage and of 
being inherently low-frequency such that they are sensitive 
to the entire moment release even for 8.5 ruptures. 
Moreover, properties such as directivity and rupture extent 
can be estimated independent of any smoothing constraints 
or other a priori infonnation [McGuire et al, 2001]. We 
quantify the extent and directivity of large ruptures using 
the second moments of a scalar source-space-time function 
describing the moment release distribution [Backus and 
Mulcahy, 1976a, 1976b; Backus, 1977a, 1977b; McGuire 
et al, 2001]. The second moments describe the length, 
width and duration of the area of greatest seismic moment 
release during an earthquake. They are defined as 

A'"’"’ = j j /(r.0(r - r(,)(r - r„fdVdt (1) 

= j J f{r,t){t-tofdVdt (2) 

A^‘’‘^ = J j fP,t){r-ro){t-to)dVdt (3) 

where / is the source-space-time function and is propor¬ 
tional to slip rate at a point on the fault, ro and fo are the 
centroid location and time, is the second spatial 

moment, is the second temporal moment, and is 
the mixed moment which describes overall rupture 
directivity [McGuire et al, 2001]. 

[lo] The second moments represent weighted averages of 
seismic moment release, and they define characteristic 
rupture dimensions that are somewhat smaller than the total 
rupture dimensions. These characteristic rupture dimensions 
are [Silver, 1983; Silver and Jordan, 1983; McGuire et al, 
2001] 

Xc{A) = 2yJhT{p(^’»yMo)n (4) 


r, = (5) 

Vc=Lclrc ( 6 ) 

v„ = (7) 

where Mq is the seismic moment and x^Cn) is the 

characteristic dimension of slip in the direction of n, which 
has its maximum value of when n corresponds to the 
largest eigenvector associated with the largest eigenvalue of 
^(2,0) xjjg characteristic duration is t^, the characteristic 
velocity is Vc and the average velocity of the instantaneous 
spatial centroid is Vq. The velocity can range from 1 to 
2 times the actual rupture velocity [McGuire et al, 2002]. A 
directivity ratio is defined by |voj/vc, such that ruptures with 
directivity ratios <0.5 are predominantly bilateral, while 
those greater than ~0.5 are predominantly unilateral. 

[11] Figure 1 demonstrates the relationship between the 
characteristic rupture dimensions and the slip distribution of 
the 2000 western Tottori earthquake. The second moments 
were calculated directly from the slip model of Iwata et al. 
[2000]. The ellipse represents the area on the fault plane that 
released the majority of the seismic moment during the 
earthquake as measured by the second spatial moment. The 
arrow indicates the rupture direetivity, which is described by 
the mixed moment. 

[ 12 ] We estimate the second moments by comparing point 
source synthetic and observed seismograms. The data 
seismogram recorded at station/component p can be 
expressed as 

Sp{x,t) = J J Gy{x,x',t — t')Mijf{x',t')djc'dt' ( 8 ) 

where Mij is the moment tensor, assumed constant during 
the earthquake, that describes fault orientation; / is the 
source-space-time function; and is the Green’s function. 
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Figure 4. Results for bilateral event 20030122 in Mexico, (a) Amplitude measurements (red) and fit 
from inversion (blue) at different frequency bands, (b) Centroid (red triangle) and characteristic rupture 
ellipses with major axes of length 0.5 (inner dashed ellipse), 1 (solid black ellipse) and 1.5 (outer 
dashed ellipse) plotted on a TPGA map. Directivity vector Vo is shown by the white arrow. The trench is 
the thin black line. The area of high seismic moment release is largely contained in a negative TPGA 
region that corresponds with the Manzanillo fore-arc basin, (c) (top) TPGA values in a single profile 
along strike of the characteristic rupture ellipse. Rupture propagates from the centroid (solid black line) 
out to both the left and right. Dashed lines mark the extent of the 1 rupture ellipse. The ends of the plot 
mark the extent of the 1.5 rupture ellipse, (bottom) Average TPGA measured over rupture ellipses of 
varying Tc. TPGA is minimized near the centroid (shown by the 0.5 Lc rupture ellipse). Rupture extent (1 
Lc) corresponds with increasing TPGA. 


The Green’s function can be expanded in a Taylor series 
about a point (x'o, t'o): 


Gfj{x,x'j 


0 


: 1 + (x' - Xj,) • V,, -h {/ - Q — 

■Gfj{x,xlt-t'o) (9) 


At low frequencies, the Taylor series can be truncated after 
the second-order terms allowing the seismogram to be 
written in terms of the zeroth, first, and second moments: 


Sp{x, t) = ^(“•“)c?^.(x, x|,, t - to)Mij -I- 
■V,Gf(x,x;,t-t')Mij 


5 of 31 

22 


















B09301 


LLENOS AND MCGUIRE: FORE-ARC STRUCTURE AND GREAT EARTHQUAKES 


B09301 


(a) 

1.5 

1 

jzi 

0.5 

0 

0 90 180 270 360 

1.5 I-'-■-^- 

0.5 |r 

3-5 mHz 

o'-^^-■- 

0 90 180 270 360 

Azimuth 






=• 


3-5 mHz 




1-3 mHz 




1.5 

1 

0.5 

0 

0 90 180 270 360 

Azimuth 


«!*!• ■ •% 


2-4 mHz 


(C) 
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Figure 5. Results for unilateral event 20060717 in Java, (a) Amplitude measurements (red) and model 
(blue), (b) See Figure 4 for symbol explanation. The area of high seismic moment release shown by the 
characteristic rupture ellipse (solid black line) stops at a positive TPGA area, (c) (top) TPGA values along 
strike of the characteristic rupture ellipse. Rupture propagates from left to right. Solid black line denotes 
the centroid location; dashed lines mark the extent of the 1 rupture ellipse, which occurs as the TPGA 
becomes more positive. The ends of the plot mark the extent of the 1.5 rupture ellipse, (bottom) 
Average TPGA measured over rupture ellipses of varying L^. TPGA is a local minimum at the centroid 
location (0.5 L^). 


+ ^ ^ - 4)Mii 

+ V,V.C?^.(x,x;,t-t')Mij (10) 

The first term of the equation on the right-hand side is the 
point source synthetic seismogram, Sp = /i*®’®*G|(x, x'o, 
t — /o)Mij- Thus the observed seismogram is described by 
the contribution from the best fitting point source perturbed 
by finite source effects. 


[ 13 ] The perturbations to the synthetic seismogram that 
describe the finite source can be estimated from the data 
using (10). To measure these anomalies, we rewrite (10) as 
Sp{\, t) = miApi{t), where m is a vector containing the 
independent elements of the zeroth, first, and second 
moments, and Api{t) is the partial derivative of the Green’s 
function specific to the ith element of m. Then the cross¬ 
correlation function of s and Sp can be expressed as 

Css{r) = Sp{t) 0Sp(t - r) fa Sp{t - t) 0 miApi(t) 

i 

= '^miSp{t - t) ® Api{t) ( 11 ) 
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Figure 6. versus duration for events in the data set, 
with a least squares fit showing the scaling between moment 
and duration. The 1994 Nicaragua and 2006 Java tsunami 
earthquakes show anomalously long durations for events of 
their magnitudes. 

A system of equations Apjiiij = bp can then be defined, 
where the measurement bp is the peak of the cross- 
correlogram Css at station p, and Apj is the cross correlation 
of Apj{t) with synthetic seismogram Sp{t — Tp). We weight 
both bp and Apj by the peak of the autocorrelogram of Sp to 
account for arbitrary differences in amplitude between 
stations. Thus, in our measurement scheme b < \ implies 
that the observed seismogram has a lower amplitude than 
the point source synthetic, as would be expected if the 
arrival is sensitive to the destructive interference associated 
with finite source effects in the frequency band being 
considered. This measurement scheme provides a straight¬ 
forward way to identify the directivity of an earthquake. 
Figure 2 shows the amplitude measurements made from a 
synthetic line source test simulating rupture propagating 
toward the east. Low amplitude measurements (6 < 1) are 
found at azimuths of around 270, away from the direction of 
propagation. 

[ 14 ] We measured frequency-dependent amplitude 
anomalies using fundamental mode Rayleigh wave data 
obtained from the GSN, GEOFON, GEOSCOPE, MEDNET, 
and China Digital seismic networks through the IRIS Data 
Management Center. The point source synthetics for an event 
were generated using the Global centroid moment tensor 
(CMT) solution for that event, available at http://www. 
globalcmt.org. We utilized two sets of point source syn¬ 
thetics: normal mode synthetics, using the one-dimensional 
(1-D) PREM Earth model [Dziewonski and Anderson, 
1981] with phase velocity maps correcting for 3-D structure 
[Ekstrom et al, 1997]; and 3-D synthetics, calculated using 
the spectral element method of Komatitsch and Tromp [1999] 
with the 3-D velocity model CRUST2.0 [Bassin et al., 2000] 
and mantle velocity model S20RTS [Ritsema and van Heijst, 
2000]. We used the 1-D normal mode synthetic seismograms 
to calculate Apj(t), the partial derivatives of the Green’s 
functions. These derivatives depend primarily on source- 
station geometry and so the accuracy of the velocity model is 


not as important as it is for the amplitude measurements, 
which we therefore made using the 3-D synthetic seismo¬ 
grams. The improvement in Earth model allowed much more 
accurate amplitude predictions to be made for the Rayleigh 
wave, especially at the higher-frequency bands used in this 
study (see Appendix A for further discussion regarding the 
use of 1-D versus 3-D synthetics). 

[15] For each event, we make measurements at a set of 
global stations with good azimuthal coverage and at a 
number of frequency bands ranging from 1 to 10 mHz. 
Waveforms are windowed around the peak of the Rayleigh 
wave using frequency-dependent window lengths prior to 
cross correlation. Stations with correlation coefficients less 
than 0.9 in any frequency band are discarded from all bands 
in the inversion to avoid errors from unmodeled heteroge¬ 
neity. The magnitude of the amplitude anomaly increases 
with frequency, however, we can only use frequency bands 
where the amplitude reduction due to finite source effects is 
less than about 60% of the point source amplitude. Above 
this band, higher-order terms in (9) become important 
(Figure 2). The useable frequency range depends on the 
spatial extent and directivity of the rupture and hence is 
different for different sized earthquakes. 

[ 16 ] The inverse problem for the second moments is 
nonunique but can be stabilized by incorporating the con¬ 
straint that the 4-D source region must have a nonnegative 
volume [Das and Kostrov, 1997; McGuire et al, 2001]. 
Other constraints are used to limit changes to the centroid 
depth as well as ensure that rupture does not occur above 
the Earth’s surface. These nonlinear constraints are 
expressed as linear matrix inequalities in the semidefmite 
programming approach of Vandenberghe and Boyd [1996]. 
The least squares objective function (11) is minimized 
subject to the various inequality constraints. The solution 
m is a 15 component vector containing the second moments 
as well as changes to the centroid time, location and seismic 
moment. 

[17] We use a leave-one-out jackknife technique to esti¬ 
mate the error of the solution [Tukey, 1984]. We divide the 
data into N subsets, where N is the number of stations. For 
the fth subset, the measurements at station i in all frequency 
bands are left out, and the inversion is performed with the 
remaining data to produce an estimate of the model param¬ 
eters. The N estimates of the model parameters then provide 
a conservative estimate of the variance of the model 
parameters [Efron and Stein, 1981]. The variance of the 
model parameters can then be used in standard error 
propagation equations to estimate the variance in derived 
quantities such as Lc and To [Bevington and Robinson, 
1992]. 

[is] The characteristic rupture ellipses for each event are 
compared to the maps of trench- parallel gravity anomaly 
(TPGA) produced by Song and Simons [2003]. To construct 
these maps, an average trench normal gravity profile for 
each subduction zone is removed from the free air gravity 
data [Sandwell and Smith, 1997] along the length of the 
subduction zone. The resulting TPGA illuminates shorter- 
wavelength features such as fore-arc basins. We compare 
the characteristic rupture dimensions with the spatial varia¬ 
tions in TPGA to test the hypothesis that a correlation exists 
between variations in the TPGA field and earthquake 
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Figure 7. Results for event 19950730 in Chile, (a) Amplitude measurements (red) and model (blue), 
(b) See Figure 4 for symbol explanation. Although the majority of the seismic moment released in this 
event occurred in a high TPGA region, the centroid is located in a local TPGA minimum, (c) (top) TPGA 
values along strike of the characteristic rupture ellipse. Rupture propagates from the centroid (black line) 
out to the limits of the rupture ellipse (dashed line). The ends of the plot mark the extent of the 1.5 Lc 
rupture ellipse. TPGA values at along-strike distances of less than 50 km should be ignored because they 
occur inland, where the TPGA measurements are not as accurate, (bottom) Average TPGA measured over 
rupture ellipses of varying L^. Inland TPGA values were masked out in calculating the average. Again a 
local minimum occurs near the centroid. 


rupture characteristics such as centroid location, rupture 
extent and directivity. 

3. Results 

[19] We estimated the second moments for 15 shallow 
thrust earthquakes (M^ > 7.5) that occurred on the plate 
interface in circum-Pacific subduction zones from 1992 to 
2006 (Figure 3). Because of the complex nature of the 
moment release during the 2004 Sumatra earthquake, this 
event was not included in our analysis. Table 1 summarizes 
the characteristic rupture dimensions of these events and 
Table 2 summarizes the second moment inversion results. 
Several representative events are discussed in greater detail 
below. Appendix B contains our Rayleigh wave measure¬ 
ments and comparisons with the TPGA field for each event. 


3.1. The 2003 Colima Mexico Earthquake 

[ 20 ] On 22 January 2003, a 7.5 earthquake occurred 
near the state of Colima, Mexieo. This event primarily 
ruptured bilaterally in the updip-downdip directions [Yagi 
et al, 2004]. Our second moment estimates (Table 1 and 
Figure 4a) support this conclusion. The inversion resulted in 
a characteristic rupture length Lc of 92 km, a duration Tc of 
29 s, and a directivity ratio of 0.44. This agrees well with 
the finite source model determined by Yagi et al. [2004] 
using a joint inversion of teleseismic body wave and strong 
motion data. 

[ 21 ] The characteristic rupture ellipse representing the 
area of greatest moment release is compared to the TPGA 
field in Figure 4b. Although the measurements and direc¬ 
tivity clearly indicate a downdip rupture, the ellipse is 
oriented more along strike. The centroid is located in a 
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Figure 8. TPGA measured along strike of rupture ellipses for each event in the data set. The solid line 
marks the centroid location, and the dashed lines mark the extent of the best fitting (1 L^) rupture ellipse. 
Rupture propagates from left to right for unilateral events and from the centroid outward in both 
directions for bilateral events. In general, the rupture ellipse limits correspond with positive changes in 
TPGA. The profiles of 19950914 and 19950730 are not entirely accurate, as part of each profile occurs 
on land, where accurate TPGA measurements are unavailable. 


highly negative TPGA region which corresponds to the 
Manzanillo basin [Wells et al, 2003]. This basin also 
overlies part of the rupture area of the 1995 8.0 Jalisco 

event (Figure B7b in Appendix B). The along-strike limits 
of the rupture ellipse are characterized by local TPGA 
maxima (Figure 4c). These results suggest that the main 
moment release filled the area underlying the Manzanillo 
basin. Figure 4c (bottom) shows the average TPGA over 
rupture ellipses of varying (with a constant aspect ratio), 
which illustrates how the TPGA changes within the rupture 
zone. The average TPGA is minimized near the centroid 
and increases as the region near the boundary of the actual 
rupture area is included (i.e., 1 in Figure 4c). Therefore 
the limits of the significant moment release are character¬ 
ized by increasing TPGA. 

3.2. The 2006 Java Earthquake 

[ 22 ] On 17 July 2006, a M^ 7.7 earthquake occurred in 
the fore arc of the Java trench. This earthquake generated a 
tsunami with a wave height of 1.8 m that struck the coast of 
Java, killing over 400 people. Preliminary results reported 


by the U.S. Geological Survey (USGS) (http://earthquake. 
usgs.gOv/eqcenter/eqinthenews/2006/usqgaf) suggest that 
this event had a similar source mechanism as that of the 
M„ 7.8 shallow thrust earthquake (also in this study) that 
occurred off the coast of Java on 2 June 1994 which also 
generated a tsunami. 

[ 23 ] For this event, we generated point source synthetics 
using the focal mechanism reported in the Global CMT 
catalog but the centroid location reported by the USGS. In 
their study of the 1994 Java earthquake, Abercrombie et al. 
[2001] found that the Global CMT locations for earthquakes 
in the Java trench appear to be systematically biased to the 
south and suggest that the 3-D velocity models used in the 
CMT inversion are inadequate in this part of the globe. 

[ 24 ] There is a very strong directivity signal at azimuths 
of around 300 (Figure 5a), indicating that rupture propagated 
unilaterally to the east-southeast. By the 3-5 mHz band, 
the amplitude of the data has dropped to 0.5 of the 
amplitude of the synthetic. Such a strong signal in this 
low of a frequency band is more typical of M^,, 8.5 events 
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Figure 9. Average TPGA over rupture ellipses of varying for each event of the data set. In general, 
the average TPGA increases with L„ indicating the centroid (reflected by the rupture ellipse of 0.5 L^) 
was located in a local TPGA minimum. The limits of the rupture ellipses which occur at 1 generally 
coincide with increases or maxima in the TPGA. 


than M.,^, 7.7 events (compare with measurements for a 
M„ 7.5 earthquake shown in Figure 4a and measurements 
for a 8.7 shown in Figure B14). The inversion for the 
second moments resulted in a length of 108 km, a 
duration of 55 s, and an average instantaneous centroid 
velocity |vo| of 1.9 km/s. The characteristic rupture length, 
relatively long duration and directivity are consistent with 
other slip models for this event [Ammon et al, 2006; Fujii 
and Satake, 2006]. The duration is surprisingly long for an 
event of this magnitude; the duration for the 7.5 Mexico 
earthquake was only 29 s. The duration and moment for 
each event in this study are shown in Figure 6, fit with a line 
illustrating the scaling of duration as the cube root of 
moment. Both this event and the 1992 Nicaragua event 
have anomalously long durations considering their magni¬ 
tudes. This shows these earthquakes belong to a class of 
slow tsunami earthquakes that radiate a large amount of 
low-frequency energy relative to their high-frequency radi¬ 
ation [e.g., Polet and Kanamori, 2000]. However, there is 
nothing obviously different about the TPGA profiles for 
these two events compared to the other events (e.g.. 
Figure 4). Thus we find no evidence in this study to suggest 


a direct relationship between structure in the TPGA field 
and tsunamigenic potential. 

[ 25 ] Figure 5b compares the characteristic rupture ellipse 
to the TPGA field. The centroid is again located in a 
negative TPGA region, and the limits of the rupture ellipse 
occur as the TPGA becomes increasingly positive 
(Figure 5c). When we consider the average TPGA over 
varying rupture dimensions, we find that the TPGA 
increases slightly near the edges of the rupture but is 
essentially flat. Hence the unilateral mpture propagation to 
the east appears to have been stopped near a local TPGA 
maximum (Figures 5b and 5c). 

3.3. The 1995 Chile Earthquake 

[ 26 ] A 8.0 earthquake occurred off the northern coast 
of Chile near the town of Antofagasta on 30 July 1995. This 
event has been the subject of a number of studies, which 
generally found that rupture initiated just south of the 
Mejillones peninsula and propagated unilaterally to the 
southwest [e.g., Deloiiis et al., 1997; Ihmle and Ruegg, 
1997; Carlo et al, 1999; Pritchard et al., 2006]. Our 
measurements and model fit are shown in Figure 7a. The 
inversion resulted in an Lc of 121 km, a Tc of 33 s, and a 
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Figure 10. (a) Rupture ellipses for events 19941228 and 20030925 off the coast of northeastern Japan 
plotted on TPGA. See Figure 4 for symbol explanation, (b) Along-strike TPGA values for (top) 
19941228 and (bottom) 20030925. Both events ruptured unilaterally, so rupture progresses from left to 
right. The 1 limit in both cases correspond with a positive TPGA slope that can be qualitatively 
correlated with transitions from velocity-weakening to velocity-strengthening behavior shown in maps 
from Miyazaki et al. [2004b] and Nishimura et al. [2004]. 


directivity ratio of 0.33, revealing a somewhat surprising 
bilateral component to the rupture. Our centroid location 
and rupture length are consistent with other slip models 
[Ihmle and Ruegg, 1997; Pritchard et al, 2006], although 
the orientation of the moment release is slightly different. 
Where other slip models are oriented predominantly along 
strike, our rupture ellipse is oriented at around a 45° angle to 
the trench. However, had our rupture ellipse been oriented 
more along strike, the overall shape of the TPGA profile 
would not change significantly, as the centroid is sur¬ 
rounded by higher TPGA regions along strike as well as 
updip and downdip (Figure 7b). 

[27] The main moment release of this event, unlike most 
of the others in this study, was located in a positive TPGA 
region, although the centroid was located in a local TPGA 
minimum (Figure 7b). These higher TPGA regions are 
associated with the Mejillones Peninsula to the north and 
the Antofagasta Ridge to the west, both of which are 
thought to be tilted fault blocks [von Huene and Ranero, 
2003]. The local TPGA minimum in which the centroid is 
located corresponds with a sediment-ponded basin behind 
the Antofagasta Ridge. The TPGA profiles show that the 
1 Lc limit occurs near the high positive TPGA region 
denoting the Antofagasta Ridge (Figure 7c). However, had 
rupture continued beyond 1 L^, the TPGA would have 
decreased and so, unlike most of the other events in the 
data set, there is no obvious correlation between positive 
changes in the TPGA field and the termination of rupture. 


3.4. Evaluation of Bias From Unmodeled Propagation 
Effects 

[ 28 ] If inadequacies of the 3-D Earth models result in a 
significant bias in our estimates of rupture area, this bias 
should be shared by earthquakes within a given subduction 
zone owing to the very long wavelength nature of the 
R1 waves we utilize. Our data set consists of events from 
essentially four regions: the northwest Pacific (six events), 
the Middle American trench (four events), the Java trench 
(three events), and the Peru-Chile trench (two events). We 
have utilized two approaches to check for bias in our 
estimated rupture areas. First, many events have been 
compared to slip distributions derived from local geodetic 
data. For instance, our rupture areas agree well with geodetic 
inversions for both the 2003 Tokachi-oki [Miyazaki et al, 
2004a] and 1995 Jalisco events [McGuire et al, 2001]. 
Moreover, there is considerable variability in the rupture 
directivity inferred for the various northwest Pacific earth¬ 
quakes indicating that systematic biases are relatively unim¬ 
portant for the majority of our events. As a further check to 
ensure that propagation effects do not bias the results, we 
compared our results to those obtained using an empirical 
Green’s function (EGF) instead of point source synthetics to 
make amplitude measurements. In general, theoretical 
Green’s functions (TGF) are preferable because of their 
superior signal-to-noise ratio at low frequencies and their 
true point source nature in time, but EGFs provide a useful 
way to check for consistency. Examples of EGF measure- 
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Figure 11. (a) Rupture ellipse (solid black ellipse) for the 2000 Tottori event plotted on the slip model 
of Iwata et al. [2000] from which it was determined. Ellipses of length 0.5 and 1.5 are also shown 
(inner and outer dashed ellipses, respectively). Centroid location is marked by the red triangle. Directivity 
is indicated by the arrow. Hypocenter location marked by the white star, (b) Rupture ellipses and 
directivity plotted on top of fracture energy Gc calculated by Dalguer et al. [2002]. (c) Gc along the axis 
of the rupture ellipse. Maxima occur near the rupture ellipse limits (dashed lines). A minimum occurs at 
the centroid (solid black line), (d) Average Gc over ellipses of varying Lc. Maximum occurs near 1 Lc- 
Centroid occurs at a local minimum. 


ments and inversion results for a northwest Pacific event 
and the two Peru-Chile trench events can be found in 
Appendix B (Figures B4, B5, and Bll). The two Peru- 
Chile trench events have ~15° and 55° rotations of their 
rupture ellipses when using EGFs versus TGFs, indicating 
that unmodeled effects may bias our estimates of rupture 
area for this subduction zone. We are not confident that the 
EGF results are better because of the large amount of scatter 
(e.g., poor SNR) at frequencies below 5 mHz for these 
events (Figure B5a and Blla). The TPGA profile changes 
qualitatively for only the 2001 Peru event. 

3.5. Summary of Results 

[ 29 ] Our survey of 15 earthquakes indicates that TPGA 
increases as a rupture approaches the edge of its eventual 
rupture patch. Figures 8 and 9 summarize the along-strike 
TPGA and average TPGA profiles for all of the events. In 
11 of the 15 events, the centroid was clearly located where 
the TPGA is a local minimum. In most cases this value was 
negative (e.g., events 19941228, 20030122), but in some 
cases it was positive (e.g., events 20050328, 19950730). 
The magnitude of TPGA variation over the rupture area also 
differs from event to event. In some cases, the average 
TPGA varied over a range of ~2 mGal, while for others 
variations occurred over a range of ~40 mGal (see events 
20060717 and 19941228 in Figure 9). This suggests that the 


relative rather than absolute values of TPGA may reflect 
rupture behavior more directly. 

[ 30 ] Our results also show that in 11 of the 15 events, the 
limits of the rupture ellipses correspond with increasing 
TPGA (Figure 9). In 3 events, the TPGA is essentially fiat 
between the centroid and the rupture edges. In one event 
(19951009), the limits correspond with decreasing TPGA, 
and so the reason for rupture to stop is not apparent in the 
TPGA field. If we exclude the South American events due 
to uncertainties in propagation effects, 9 of 13 events have 
increasing TPGA. Thus, in the majority of the events in our 
data set, the 1 Lc limit corresponds with a positive gradient 
in TPGA (Figure 9). This suggests that a relationship exists 
between the TPGA field and physical conditions on the 
plate interface that control the extent of individual ruptures. 

4. Discussion 

[ 31 ] In agreement with previous studies [Song and 
Simons, 2003; Wells et al, 2003], our results demonstrate 
that the areas of greatest seismic moment release for large 
subduction zone earthquakes tend to occur beneath local 
minima in the TPGA field. Typically these minima denote 
the presence of a fore-arc sedimentary basin. Additionally, 
we demonstrated that the limits of the two-dimensional 
rupture areas correlate with positive gradients of TPGA 
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Figure Al. (top left) Map showing great circle paths from the source of the 2003 Mexico earthquake to 
stations. Red stations were used in the final inversion. Comparison of 1-D and 3-D synthetic waveforms 
at stations MHV, LCO, and ATD (blue triangles on the map) to actual data in the 10-12 mHz band. The 
amplitude measurements are markedly different, especially at LCO. 


(Figure 9). Our observations lead to two primary conclu¬ 
sions regarding the connection between TPGA and the 
frictional conditions on the plate interface. First, the relative 
rather than absolute values of TPGA have the greatest 
correlation to individual rupture behavior. Second, there 
appears to be a correlation between the resistance to rupture 
growth (either spatial variations in current shear stress level 
or material properties) and the edges of local TPGA 
minima. Although these minima often correlate with fore¬ 
arc basins, they do not identify them uniquely because of 
other variations in seafloor topography and density which 
contribute to the gravity field. However, Wells et al. [2003] 
demonstrated that the basin edges are spatially correlated 
with maxima in gravity gradient rather than absolute values 
of the free air gravity field, which they then correlate with 
updip and downdip frictional stability transitions. In this 
section, we further explore the potential mechanisms con¬ 
necting the upper plate fore-arc basins reflected in the 
TPGA field and the resistance to rupture propagation 
observed to occur near the edges of these basins. 

4.1. Fore-Arc Basin Formation and Wedge Mechanics 

[ 32 ] A number of 2-D models have examined the link 
between fore-arc basin formation and conditions along the 


plate interface. Byrne et al. [1988] suggested that the 
presence of fore-arc basins indicates strong material at depth 
which causes a change in the critical taper angle of the 
wedge and allows backstops and outer-arc highs to fonn 
which trap sediments. Recent numerical work has explored 
the mechanics of the upper plate that allow the formation of 
fore-arc basins {Fuller et al., 2006; Wang and Hu, 2006]. 
Wang and Hu [2006] divide the fore arc into an undeform¬ 
ing inner wedge and an actively deforming outer wedge. 
The inner wedge provides a stable platform for fore-arc 
basins to form on and overlies the velocity-weakening 
seismogenic zone, while the outer wedge overlies the 
velocity-strengthening region updip of the seismogenic 
zone. Fuller et al. [2006] also suggest that fore-arc basins 
tend to form in stable wedges, and the lack of defonnation 
in such wedges could mean that permeability is reduced in 
the direction normal to the plate interface, allowing process¬ 
es such as thermal pressurization to generate elevated pore 
pressures that could significantly weaken the fault [Wibber- 
ley and Shimamoto, 2005]. Song and Simons [2003] infer 
that strongly negative TPGA values correlate with increases 
in the shear traction on the plate interface. Large shear 
tractions tend to result in more stick-slip behavior [Marone, 
1998]. While the existing models are fairly simplistic, they 
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Figure A2. Amplitude measurements in different frequency bands made using 1-D synthetics (blue) 
and 3-D synthetics (red). At higher frequencies, the measurements become significantly different, 
particularly in northern azimuths. There is also more scatter in the 1-D synthetic measurements. 


demonstrate mechanical reasons to expect fault strength 
and/or fluid pressure to be different under fore-arc basins 
than underneath adjacent areas. In general, these frictional 
differences do not appear to be related to variations in the 
thermal structure of the plate interface. 

4.2. How Great Earthquakes Stop: Dynamic Versus 
Quenched Heterogeneity 

[33] The resistance to rupture growth that causes earth¬ 
quakes to stop can arise from along-strike variations in 
shear stress (dynamic heterogeneity) or along-strike varia¬ 
tions in fault frictional properties (quenched heterogeneity) 
[Guatteri and Spudich, 2000; Shaw, 2000]. In numerical 
models, time-dependent, dynamic heterogeneities can pro¬ 
duce complex earthquake sequences even on faults with 
uniform frictional properties owing to the low-stress regions 
leftover from past ruptures [Shaw, 2000]. However, the 
locations of these dynamic heterogeneities vary on the scale 
of an earthquake cycle (100 s of years). In contrast, the fore¬ 
arc basins (and associated TPGA anomalies) reflect the 
average stress state on the plate interface and within 
the wedge over timescales of basin formation (at least a 
million years) [Fuller et al, 2006]. Thus, if large ruptures 
are stopped simply by encountering the low-stress regions 
leftover from previous ruptures, as opposed to some long- 
lived material heterogeneity, then averaging over many 
earthquake cycles (or equivalently many subduction zones) 
should remove any correlation with the TPGA field. Our 
observed correlation between TPGA variations and the 
boundaries of individual ruptures requires that long-lived, 
upper plate structural/geological heterogeneity is a first- 
order control on along-strike rupture extent. 

[34] One type of along-strike frictional heterogeneity on a 
fault that may prevent large throughgoing ruptures are 
patches of velocity-strengthening material that exist due to 


compositional or thermal anomalies. The correlation 
between positive TPGA gradients and a transition from 
velocity-weakening to velocity-strengthening behavior may 
be qualitatively assessed in a region such as northeastern 
Japan. Owing to the dense instrumentation network, a 
number of studies have mapped the long-term seismogenic 
behavior of the plate interface in this region [Miyazaki et al., 
2004b; Nishimura et al, 2004; Yamanaka and Kikuchi, 
2004]. Two earthquakes in our study occurred here: the 1994 
Mw 7.7 Sanriku-oki earthquake off the coast of northern 
Honshu (19941228) and the 2003 M.,^, 8.3 Tokachi-oki 
earthquake off the coast of Hokkaido (20030925). In their 
study of the 2003 earthquake, Miyazaki et al. [2004b] 
inferred that the afterslip of this event indicated velocity¬ 
strengthening behavior and found that it occurred in areas 
surrounding the greatest moment release (both along strike 
and updip). Significant afterslip also occurred in the areas 
surrounding the greatest coseismic slip of the 1994 Sanriku- 
oki earthquake [Nishimura et al., 2004]. 

[35] Our characteristic rupture ellipses for these two 
events are shown in Figure 10a. For 19941228, the limits 
of the rupture ellipse occur when the TPGA becomes 
positive. The positive TPGA regions correspond relatively 
well with the weak seismic coupling region shown by 
Nishimura et al. [2004]. This observation seems to suggest 
a correlation between positive TPGA values and velocity¬ 
strengthening behavior. However, the rupture ellipse for 
20030925 shows that the rupture failed to fill the entire 
basin and stopped in a negative TPGA region. This suggests 
the presence of velocity-strengthening materials even 
though the TPGA remains negative. Instead, the transition 
from velocity-weakening to velocity-strengthening behavior 
seems to be correlated with positive TPGA gradients. 
Beyond the limits of the rupture ellipses for both events, 
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Figure Bl. Measurements and results for unilateral event 19940602 in Java. See Figure 4 for symbol 
explanation. 


the TPGA increases (Figure 10b). This change in the TPGA 
corresponds well in location with changes in the seismic 
coupling seen in maps drawn by Nishimum et al. [2004] and 
Miyazaki et al. [2004b] for the 1994 and 2003 earthquakes, 
respectively. At least in the Tokachi-oki region, there appears 
to be a link between increases in TPGA and transitions 
from velocity-weakening to velocity-strengthening behavior. 

[ 36 ] A second type of frictional heterogeneity that can 
play a role in stopping earthquake rupture is regions of 
velocity weakening material with very high fracture energy. 
In subduction zones, this may represent areas of increased 
fault roughness due to the bathymetry of the incoming plate. 
Figure lib shows the characteristic rupture ellipse for the 
2000 Tottori earthquake plotted on fracture energy calcu¬ 
lated by Dalguer et al. [2002]. High fracture energies are 
associated with the limits of the rupture ellipse (Figures 11c 
and lid) suggesting that the region near 1 identifies the 
location of significant resistance to rupture propagation that 


was necessary to arrest a particular rupture. The fracture 
energy profiles in Figures 11c-lid are remarkably similar 
to the typical TPGA profiles found in this study (Figures 8-9) 
despite being nominally unrelated variables. Thus the var¬ 
iations we observe in TPGA around the 1 Lc portion of the 
rupture zone may reflect a surprisingly direct connection 
between the edges of sedimentary hasins and the variations 
in frictional properties of the plate interface. 

[ 37 ] The primary difference between the fracture energy 
and frictional stability transitions as explanations for the 
boundaries of M„ 8 ruptures lies in the role of aseismic slip. 
In the fracture energy case, most of the plate motion would 
be expected to be made up by abundant smaller (M.,^ <8) 
earthquakes, while in the stability transition explanation, 
much of the plate motion in these regions could occur 
aseismically, likely as afterslip. Thus distinguishing 
between these two types of frozen heterogeneity will require 
geodetic observations over multiple earthquake cycles, hut 
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Figure B2. Measurements and results for unilateral event 19941228 in Japan. See Figure 4 for symbol 
explanation. 


at least for the Tokachi-oki region, along-strike variations in 
frictional stability are important [Miyazaki et al, 2004b]. 

5. Conclusion 

[ss] The spatial heterogeneity of frictional properties 
which cause variations in seismogenic behavior along a 
subduction zone can be linked to variations in fore-arc 
TPGA and geological structure [Song and Simons, 2003]. 
To further investigate this relationship, we estimated the 
rupture charaeteristics of a global data set of large 
shallow subduction zone earthquakes and compared them 
with the TPGA field. The centroid locations typically 
correspond with local TPGA minima, and the limits of 


the areas of significant moment release correspond with 
positive TPGA gradients (Figure 9). These gradients 
appear to be linked to changes in the frictional properties 
along the plate interface that control rupture behavior. 
Owing to the inherently long timescales required for fore¬ 
arc basin formation, the correlation between the TPGA 
field and rupture termination regions indicates that long- 
lived material heterogeneity rather than short timescale 
stress heterogeneities are responsible for arresting most 
great subduction zone ruptures. Variations in TPGA may 
therefore be used to not only detennine the long-term 
seismogenic behavior along the strike of subduction zones 
but also to provide estimates of the boundaries of 
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Figure B3. Measurements and results for unilateral event 19950914 in Mexico. See Figure 4 for symbol 
explanation. 


individual future ruptures for probabilistic seismic hazard 
analysis. 

Appendix A: One-Dimensional Versus Three- 
Dimensional Synthetics 

[39] In this appendix, we discuss the reasons for using 
both 1-D and 3-D synthetic seismograms in this study. We 
use 1-D nomal mode synthetic seismograms to calculate 
the partial derivatives used in the inversion and 3-D 
synthetics to make the actual amplitude measurements. 
The 1-D synthetics are much faster to compute than the 
3-D synthetics, which require a 64-processor Linux cluster, 
and so they are considerably more efficient to use for the 
numerical calculations of the partial derivatives of the 
Green’s functions. The use of 3-D synthetics rather than 
1-D synthetics to calculate the partial derivatives in (10) 
would make little difference in the inversion results. The 
partial derivatives are calculated from perturbations in 
source location, and so are more sensitive to source-station 
geometry than to the velocity model used to generate the 


synthetics. The results from a test event confirm that little is 
gained by using 3-D synthetics for the partial derivatives; 
therefore, because the 3-D synthetics require much more 
computational effort, we use the 1-D synthetics to calculate 
the partial derivatives. 

[ 40 ] However, the 3-D synthetics produce more accurate 
R1 amplitude predictions at higher frequencies, so they are 
used to make the amplitude measurements. The difference 
in amplitude predictions can be seen by comparing the 
1-D and 3-D synthetics with the seismograms observed 
at three stations at different azimuths for the 2003 M„ 7.5 
earthquake in Mexico (Figure Al). Amplitude measure¬ 
ments made using 1-D synthetics would have been lowest 
at station MHV (az = 21). With the 3-D synthetics, the 
lowest measurement occurs at station LCO (az = 145), 
which indicates a completely different direction of rupture 
directivity. 

[ 41 ] The magnitude of amplitude anomaly measured with 
1-D and 3-D synthetics can also be very different (see 
station LCO in Figure Al). Figure A2 compares the ampli¬ 
tude measurements made using 1-D and 3-D synthetics for 
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Figure B4. Measurements and results for event 19941228 in Japan made using a 6.7 aftershock that 
occurred on 6 January 1995 as an empirical Green’s function, (a) Amplitude measurements (red), fit from 
inversion (blue), and measurements from stations discarded in EGF due to poor SNR but retained in the 
original analysis (black). Note scatter in the lowest bandpass, as well as increase in amplitude with 
frequency due to differences in centroid depth of the EGF and the main shock, (b) See Figure 4 for 
symbol explanation. Compared to the original ellipse, this rupture ellipse has rotated by ~28°, although 
the centroid location, rupture length, and directivity are similar, (c) See Figure 4 for symbol explanation. 
Qualitatively, these TPGA plots changed very little from the original analysis. 
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Figure B5. Measurements and results for event 19950730 in Chile using a 5.9 aftershock that 
occurred on 3 August 1995 as an EGF. (a) See Figure B4 for symbol explanation. Note the large amount 
of scatter in the 2-4 mHz band, (b) See Figure 4 for symbol explanation. Compared to the TGF results, 
the ellipse has rotated by ~14° but remains oblique to the strike of the trench, (c) See Figure 4 for symbol 
explanation. Although these plots have changed slightly from the original results, the overall conclusion 
remains the same qualitatively. 
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Figure B6. Measurements and results for unilateral event 19950914 in Mexico. See Figure 4 for symbol 
explanation. 
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Figure B7. Measurements and results for unilateral event 19951009 in Mexico. See Figure 4 for symbol 
explanation. 
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Figure B8. Measurements and results for bilateral event 19951203 in the Kurile Islands. See Figure 4 
for symbol explanation. 
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Figure B9. Measurements and results for unilateral event 19971205 in Kamchatka. See Figure 4 for 
symbol explanation. 
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Figure BIO. Measurements and results for bilateral event 20010623 in Peru. See Figure 4 for symbol 
explanation. 
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Figure Bll. Measurements and results for event 20010623 in Peru, using a M„ 6.7 aftershock that 
occurred on 26 June 2001 as an EGF. (a) See Figure B4 for symbol explanation. Again, a large amount of 
scatter appears in the lowest frequency band, (b) See Figure 4 for symbol explanation. Compared to the 
TGF results, the ellipse has rotated by ~54°; the orientation of the ellipse for this event is not well 
constrained, although the centroid location is. The directivity ratio remains very low. (c) See Figure 4 for 
symbol explanation. The rotation of the ellipse has caused these plots to change significantly. 
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Figure B12. Measurements and results for unilateral event 20030925 in Japan. See Figure 4 for symbol 
explanation. 
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Figure B13. Measurements and results for unilateral event 20031117 in the Rat Islands. See Figure 4 
for symbol explanation. 
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Figure B14. Measurements and results for unilateral event 20050328 in Sumatra. See Figure 4 for 
symbol explanation. 
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Figure B15. Measurements and results for unilateral event 20061115 in the Kurile Islands. See Figure 4 
for symbol explanation. 
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complete station sets. Significant differences in the meas¬ 
urements occur at northerly azimuths in the highest fre¬ 
quency band. These differences alter the results of the 
inversion. An inversion with the 1-D synthetics resulted in 
a rupture length of 51 km. Using the 3-D synthetics in the 
inversion resulted in a rupture length of 92 km. Therefore, 
because the improvement in propagation modeling is sig¬ 
nificant enough to affect the inversion results, we use 3-D 
synthetics for the amplitude measurements. 

Appendix B: Results for Additional Events 

[ 42 ] This section contains the results for each event in our 
data set (Figures B1-B15; see Table 1 for list of events). 
For each event, we show Rayleigh wave measurements, 
estimated rupture area plotted on a TPGA map, and TPGA 
profiles. For three events (19941228, 19950730, and 
20010623), the results using measurements made with 
empirical Green’s functions as described in section 3.4 are 
also shown (Figures B4, B5, and Bll, respectively). These 
results can be directly compared with those obtained by 
using point source synthetics (theoretical Green’s functions) 
as our analysis does (Figures B3, 7, and BIO). 
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software freely distributed by Wessel and Smith [1998]. A. Llenos was 
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Chapter 3: Modeling seismic swarms triggered by aseismic transients 


Abstract 

The rate of earthquake occurrence varies by many orders of magnitude in a given 
region due to variations in the stress state of the crust. Our focus here is on variations in 
seismicity rate triggered by transient aseismic processes such as fluid flow, fault creep or 
magma intrusion. While these processes have been shown to trigger earthquakes, 
converting observed seismicity variations into estimates of stress rate variations has been 
challenging. Essentially aftershock sequences often obscure changes in the background 
seismicity rate resulting from aseismic processes. Two common approaches for 
estimating the time dependence of the underlying driving mechanisms are the stochastic 
Epidemic Type Aftershock Sequence model (ETAS) [Ogata, Y., (1988), Statistical 
models for earthquake occurrences and residual analysis for point processes, J. Am. Stat. 
Assoc., 83, 9-27.] and a physical approach based on the rate- and state-model of fault 
friction [Dieterich, J., (1994), A constitutive law for rate of earthquake production and its 
application to earthquake clustering, J. Geophys. Res., 99, 2601-2618.]. The models have 
different strengths that could be combined to allow more quantitative studies of 
earthquake triggering. To accomplish this, we identify the parameters that relate to one 
another in the two models and examine their dependence on stressing rate. A particular 
conflict arises because the rate-state model predicts that aftershock productivity scales 
with stressing rate while the ETAS model assumes that it is time independent. To resolve 
this issue, we estimate triggering parameters for 4 earthquake swarms contemporaneous 
with geodetically observed deformation transients in various tectonic environments.We 
find that stressing rate transients increase the background seismicity rate without 
affecting aftershock productivity. We then specify a combined model for seismicity rate 
variations that will allow future studies to invert seismicity catalogs for variations in 
aseismic stressing rates. 


Published as A. L. Llenos, J. J. McGuire, and Y. Ogata, Modeling seismic swarms triggered by 
aseismic transients. Earth Planet. Sci. Lett., 281, 59-69, doi: 10.1016/j.epsl.2009.02.011, 2009. 
Reprinted with permission of Elsevier, 2009. 
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The rate of earthquake occurrence varies by many orders of magnitude in a given region due to variations in 
the stress state of the crust. Our focus here is on variations in seismicity rate triggered by transient aseismic 
processes such as fluid flow, fault creep or magma intrusion. While these processes have been shown to 
trigger earthquakes, converting observed seismicity variations into estimates of stress rate variations has 
been challenging. Essentially aftershock sequences often obscure changes in the background seismicity rate 
resulting from aseismic processes. Two common approaches for estimating the time dependence of the 
underlying driving mechanisms are the stochastic Epidemic Type Aftershock Sequence model (ETAS) 
[Ogata, Y., (1988), Statistical models for earthquake occurrences and residual analysis for point processes, 
J. Am. Stat. Assoc., 83, 9-27.] and a physical approach based on the rate- and state-model of fault friction 
[Dieterich, J., (1994), A constitutive law for rate of earthquake production and its application to earthquake 
clustering,]. Geophys. Res., 99, 2601-2618.]. The models have different strengths that could be combined to 
allow more quantitative studies of earthquake triggering. To accomplish this, we identify the parameters that 
relate to one another in the two models and examine their dependence on stressing rate. A particular conflict 
arises because the rate-state model predicts that aftershock productivity scales with stressing rate while the 
ETAS model assumes that it is time independent. To resolve this issue, we estimate triggering parameters for 
4 earthquake swarms contemporaneous with geodetically observed deformation transients in various 
tectonic environments. We find that stressing rate transients increase the background seismicity rate without 
affecting aftershock productivity. We then specify a combined model for seismicity rate variations that will 
allow future studies to invert seismicity catalogs for variations in aseismic stressing rates. 

© 2009 Elsevier B.V. All rights reserved. 


1. Introduction 

Earthquake swarms are time periods of elevated seismicity rate 
that lack an obvious mainshock, and they are one of the clearest 
signals that many processes in the crust have variations on time scales 
of hours to days. The most common time periods of increased 
seismicity rate are the aftershock sequences that follow all large 
crustal earthquakes and generally decay away according to Omori's 
empirical law (Utsu, 1961). The term swarm has been used 
qualitatively by seismologists for nearly a century to describe 
temporal clusters of earthquakes that are not well described by 
Omori's law (Richter, 1958). Swarms are common in volcanic regions 
and have been explained as resulting from the stress perturbations 
during magma intrusions (e.g., Einarsson and Brandsdottir, 1980; 
Dieterich et al., 2000; Smith et al., 2004) as well as from the 
movements of volatiles such as CO 2 (e.g.. Prejean et al„ 2003; Hainzl 
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and Ogata, 2005). Similarly, earthquake swarms are common in 
regions of aqueous fluid flow such as geothermal areas (Hill et al., 
1975) and during hydrofracture experiments (Audigane et al., 2002; 
Shapiro et al., 2005; Bourouis and Bernard, 2007). Thus, a clear 
intuition has developed that swarms are driven by aseismic events 
that temporarily modify the stress state within the crust. Toda et al. 
(2002) recently formalized this hypothesis for a swarm of -7000 
earthquakes in the Izu islands that was associated with a large dike 
intrusion. They demonstrated that the seismicity rate varied spatially 
in proportion to the variations in the stress rate increase caused by the 
magmatic intrusion. 

Recently a number of earthquake swarms have been found in 
association with times when a fault undergoes a large amount of slip 
without radiating seismic waves. These events are often termed slow 
earthquakes, silent earthquakes or creep events and require high 
quality geodetic data to detect owing to their lack of seismic radiation. 
Swarms triggered by such aseismic fault slip have been found in a 
number of tectonic regions. Ozawa et al. (2007) found swarms 
coincident with repeating slow earthquakes on the subduction zone 
thrust interface offshore of central Honshu. In these cases the slow 
event typically has a moment magnitude of M„~6.5 while the largest 
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earthquakes in the swarm are Mw~4, indicating that the vast majority 
of fault motion happens aseismically and only a few small patches fail 
seismically. A similar behavior was observed for a continental strike- 
slip fault by Lohman and McGuire (2007) in the Salton Trough region 
of California, where a swarm of -1000 JV/w<5.1 earthquakes was 
triggered by a 5.7 slow event. Wolfe et al. (2007) found swarms of 
-10-50 earthquakes associated with slow events on a detachment 
fault beneath Kilauea's south flank. Similarly a swarm of -1700 
earthquakes in a volcanic region in Japan was also found to be 
associated with a much larger aseismic slip transient on reverse faults 
(Takada and Furuya, in review). Collectively these studies indicate that 
relatively modest earthquake swarms with events in the magnitude 
4-5 range often result from much larger aseismic slip transients that 
generate microseismicity by loading neighboring regions of a fault 
system. Additionally, surveys of seismicity catalogs by Vidale et al. 
(2006) and Vidale and Shearer (2006) have found that swarms are 
widespread phenomena in California and japan and often cover an 
unusually large area for their cumulative seismic moment, a property 
that corresponds well with the low stress drops observed for shallow 
aseismic creep events (Brodsky and Mori, 2007). 

It is difficult to untangle the contribution of any time-dependent 
driving process from an earthquake catalog because of the preponder¬ 
ance of standard earthquake-earthquake triggering (e.g., aftershock 
sequences). Given a triggering model that utilizes an aftershock 
triggering exponent a, the Gutenberg-Richter parameter b, and an 
offset AMafter describing the magnitude difference between a mainsbock 
and its largest probable aftershock, the branching ratio n = / 

(b — a) describes the average fraction of a catalog consisting of triggered 
earthquakes (Helmstetter and Sornette, 2002; Boettcher and Jordan, 

2004) . Using values ofa w 0.8, fa 1 and AMafter« 0.9 that are consistent 
with southern California seismicity data (Helmstetter, 2003; Helmstetter 
and Sornette, 2003), up to roughly 90% of earthquakes in this region are 
triggered by other earthquakes. This number, however, is highly 
dependent on the value for a, which is likely between 0.8 and 1, and the 
parameter AMaSer can also range from 0.9 to 1.2 (Helmstetter, 2003; 
Helmstetter and Sornette, 2003). Yet even with these ranges of values, 
around 60-90% of earthquakes in the catalog are aftershocks. Thus, even 
when some aseismic process is triggering an elevated rate of seismicity, 
that seismicity will generate its own aftershock sequences, which will 
ultimately comprise a significant fraction of the earthquake catalog. There 
are currently two main approaches to analyzing seismicity rate variations; 
stochastic models such as the Epidemic Type Aftershock Sequence (ETAS) 
model (Ogata, 1988), and physical models such as the rate- and state- 
dependent friction model (Dieterich, 1994). 

The ETAS stochastic model is an effective way to detect anomalous 
seismicity rates. By modeling earthquake occurrence as a point 
process described by just a few optimizable parameters, the model 
can detect time periods that are not well described by a stationary 
stochastic process (McGuire et al., 2005). Recently, studies have 
utilized a space-time version of ETAS to relate non-stationary 
seismicity rates to regional stress changes (Ogata, 1998, 2004, 

2005) . The difficulty with this approach is that it lacks a quantitative 
relationship between seismicity rate variations and stress/stressing- 
rate variations. However, it has been used to resolve time dependence 
of the background triggering rate by binning unusually large earth¬ 
quake swarms into smaller (moving window) time periods that are 
assumed to have a stationary background rate within the time 
window (Hainzl and Ogata, 2005). 

An alternative approach that is being utilized to map seismicity 
rate variations directly into stressing rate variations is a physical 
model based on rate- and state-dependent friction (Dieterich, 1994; 
Dieterich et al., 2000). This model incorporates several properties of 
laboratory friction measurements including an Omori-like response to 
a step change in stress-level. It has had several successful applications 
including retrieving the magnitude of stress steps using laboratory 
derived friction parameters (Dieterich et al., 2000) and predicting the 


spatial distribution of seismicity rate changes and aftershock sequence 
durations based on a geodetically derived model of stressing rate 
transients (Toda et al., 2002). However, both these applications 
occurred in volcanic regions where aftershock sequences are often 
subdued due to high geothermal gradients (Kisslinger and Jones, 1991; 
Ben-Zion and Lyakhovslcy, 2006). In contrast, Toda and Matsumura 
(2006) used this method to estimate spatio-temporal stress changes 
from seismicity rate changes during a Mw 7.0 slow subduction zone 
earthquake in Tokai, Japan. Despite the extremely large magnitude of 
the slow event, the stressing rate changes retrieved by the rate-state 
inversion were not clearly distinguishable from other variations in the 
area. Some of this lack of resolution likely results from the 
contamination of moderate aftershock sequences in the stress vs. 
time curves produced by the rate-state inversion algorithm. 

We seek to combine tbe strengths of the ETAS and rate-state 
approaches to develop an effective tool to detect anomalous seismicity 
rates and relate them to changes in stressing rates caused by physical 
processes. There have been recent attempts to combine the two 
models of seismicity rate for different purposes. For example. Console 
et al. (2006,2007) combine the two models in order to produce a new 
model of earthquake clustering that incorporates physical constraints 
with a minimum number of free parameters. However, a combined 
ETAS/rate-state model that can be used in a single algorithm to detect 
anomalous stressing rates from seismicity rates has not been 
developed yet. Ogata (2005) demonstrated that even small changes 
in stress can cause anomalies in seismicity rate, and so a combined 
ETAS/rate-state model of seismicity rate has the potential to be a 
highly sensitive detector of transient deformation. 

To develop such a combined model of seismicity rate, we first 
identify parameters that are related between the two models and 
examine their dependence on stressing rate. To clarify the stressing 
rate dependence of the aftershock parameters, we analyze data from 4 
different earthquake swarms in various tectonic settings. We then 
specify a functional form for the seismicity rate in a combined ETAS/ 
rate-state model, in which aseismically-triggered and coseismically- 
triggered components of seismicity rate are independent of one 
another. This suggests that an aseismically-triggered seismicity rate 
can be isolated from a catalog and used to directly estimate stressing 
rate changes associated with transient deformation. 

2. Models 

In general, the seismicity rate R in a catalog is a function of the 
stressing rate S acting on a fault (Dieterich, 1994). Earthquake catalogs 
contain seismicity triggered by different underlying mechanisms, 
such as earthquake-earthquake interactions, aseismic deformation, or 
background plate tectonic motion. Therefore, in our model, we 
consider three primary contributions to R: 

R{x,t) =f(S} =f(R^,Rc,Ry) 

where R^ represents the seismicity rate triggered by aseismic 
processes such as slow slip or dike intrusion, Rc reflects seismicity 
triggered by other earthquakes (e.g., aftershock sequences), and Rj is 
triggered by long-term tectonic loading. Because tbe tectonic 
component Rt is presumably small if aseismic processes are occurring, 
we simply combine it with Ra, so that R(x, t)a^f(R/),,Rc). 

In order to develop a model that can quantitatively relate 
stressing rates to seismicity rates, we need to know the stressing 
rate dependence of each of the components of R. Toda et al. (2002) 
and Lohman and McGuire (2007) showed that seismicity rates 
clearly increase during periods of increased stressing rate caused by 
aseismic processes. Both studies found that the increase in earth¬ 
quake rate was approximately equal to the increase in stressing rate 
so Ra likely scales linearly with stressing rate as predicted by the 
rate-state model. 
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Fig. 1. Calculation of seismicity rates using forward simulation of rate-state model equations for two different stress histories, with Acj= 0.1 MPa. a) Stress histories for two cases: 
i) increase in stressing rate by 2 orders of magnitude (blue), and ii) increase in stressing rate by 2 orders of magnitude plus sudden stress changes (earthquakes) (red), b) Stressing 
rate histories for both cases, c) Evolution of y from using Eqs. (3) and (4). d) Seismicity rate estimates obtained from stressing rates using Eq. (2). 


The stressing rate dependence of Rc is less certain. The two main 
approaches to modeling the change in seismicity rate following an 
earthquake (i.e., Rc) are the stochastic Epidemic Type Aftershock 
Sequence (ETAS) model (Ogata, 1988), which is based on Omori's law 
(Omori, 1894), and the physically based rate-state friction model, 
which reproduces Omori's law following a sudden stress change 
(Dieterich, 1994). In this section we summarize the two models and 
compare how they predict Rq changes with stressing rate. 

2.1. ETAS model 


The ETAS model is a point process model that generalizes the 
modified Omori law (Omori, 1894; Utsu, 1961; Ogata, 1988). In this 
model, every aftershock has some probability of generating its own 
aftershocks. Therefore, the seismicity rate at some time t can be obtained 
by summing the aftershock sequences produced by each event that has 
occurred prior to time t plus a background seismicity rate p: 


R{t) — p + 


E 


(t-tj + C)P 


( 1 ) 


where c and p are the Omori decay parameters, a is related to the 
efficiency of an earthquake of a given magnitude at generating aftershocks, 
and K reflects the aftershock productivity of a mainshock. These 
parameters are generally obtained using maximum likelihood estimation 
from the observed times q and magnitudes Mi of earthquakes in a catalog, 
given the magnitude of completeness Me of the catalog (Ogata, 1988). The 
summation in Eq. (1) is essentially the coseismic component of seismicity 
rate (Rc), as it contains all the aftershock sequences in the catalog. Because 
the ETAS parameters are not explicitly related to stressing rate, Rc also is 
independent of stressing rate in the ETAS model. 


2.2. Rate- and state-dependent friction model 


state-dependent friction can be linked to a stressing rate Sthrough the 
following equations that assume normal stress is constant (Dieterich, 
1994): 



where r is the steady-state seismicity rate for the same magnitude 
interval associated with a reference stressing rate S,, S is a modified 
Coulomb stress function, 7 is a state variable, and A is a fault 
constitutive parameter. We assume that the normal stress a remains 
constant, and as a result treat Atr as a constant frictional parameter 
and S as a shear stressing rate. 

Using this formulation, given some knowledge of regional back¬ 
ground seismicity and Act, the stressing rate can be obtained from an 
observed seismicity rate simply by integrating Eqs. (2) and (3) 
(Dieterich et al., 2000). Eig. 1 illustrates the relationship between 
stressing rate and seismicity rate in the rate-state model. Given a 
stress history that involves a change in stressing rate by two orders of 
magnitude (Eig. la-b, blue), Eq. (3) can be used to calculate the 
related change in 7 (Fig. Ic, blue), which can then be used in Eq. (2) to 
obtain the change in seismicity rate. Fig. Id demonstrates that 
following a change in stressing rate by a factor of 100 , a similar 
change in seismicity rate occurs after a time lag that is related to the 
parameter Act. 

Now consider a stress history that includes the same change in 
stressing rate as well as earthquakes (sudden stress steps) (Fig. la-b, 
red). Dieterich (1994) derived a solution for 7 given a sudden stress 
change (Eq. Bll in the 1994 paper): 


7 = 70 exp 


-AS" 

Act 


In the rate-state model, the seismicity rate R for a given magnitude 
interval observed on a population of faults governed by rate- and 
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Using both Eqs. (3) and (4), the evolution of y associated with this 
stress history can also be determined (Fig. Ic, red) and Eq. (2) used to 
obtain the seismicity rate (Fig. Id, red). 

This simple case shows that the rate-state model predicts that 
certain parameters describing aftershock decay change with stressing 
rate. For example, consider two of the sudden stress steps (earth¬ 
quakes) shown in Fig. lb, the first (at time f = 100) which occurs prior 
to the stressing rate change and the third (at time f=1150) which 
occurs well after the stressing rate change. The seismicity rate 
following the first earthquake, which occurs during low stressing 
rates, takes longer to decay to the background rate than that following 
the second earthquake, which occurs during high stressing rates 
(Fig. Id); thus, this characteristic relaxation time, depends on 
stressing rate as seen in the Miyake-Jima swarm (Toda et al., 2002). 

This case also demonstrates that the change in seismicity rate 
immediately following the stress step also varies with stressing rate. 
Although the first and second earthquakes both had the same stress 
change, the peak seismicity rate is higher for the second earthquake 
than for the first. Because the change in seismicity rate AR due to the 
stress change depends on the change in y, this is a direct consequence 
of Eq. (4). If we define Ay to be the change in y due to the stress step 
(i.e.. Ay=y — yo), it is easy to see from Eq. (4) that Ay will be a function 
of yo (i.e., the value of y prior to the stress step). Since the second 
earthquake occurs during the higher stressing rate, it has a lower yo, 
and therefore a smaller Ay and a larger AR than the first earthquake 
(Fig. Ic). Thus the seismicity rate immediately following a stress step 
depends on the stressing rate. 

This prediction of the rate-state model can also be shown with a 
second simulation. Assuming a constant stressing rate following a stress 
step and steady-state seismicity rate prior to a stress step, Dieterich 
(1994) used Eqs. (2)-(4) to express the seismicity rate following a 
stress step as: 


R(t) = 




-exp I 




-1 


S 5^ 0 


exp 


-F 1 


(5) 


where 

t,=A(T/S (6) 

is the characteristic relaxation time related to the time it takes for the 
seismicity rate to return to background levels. This takes the form of 
Omori's law for t<t^. 

Eq. (5) can be used to compare the seismicity rate change due to a 
uniform stress step during different magnitudes of stressing rate. 
Because Eq. (5) is only valid when the seismicity rate prior to the stress 
step is at steady-state, we assume that the reference seismicity rate r 
(i.e., the seismicity rate prior to the stress step) has already achieved a 
steady-state value at each of the stressing rate levels. Therefore, r will 
have a different value at each stressing rate (Fig. 2a, dashed lines), 
because as Eqs. (2) and (3) demonstrate, the steady-state seismicity 
rate scales with stressing rate (Fig. 1 ). Furthermore, we assume that the 
stressing rate before and after the stress step remains constant (i.e., 
Sr = S). Given these assumptions, we can use Eq. (5) to predict the 


Fig. 2. a) Seismicity rate R calculated with Eq. (5) at different magnitudes of stressing 
rate S relative to a background stressing rate Sb = 0.1 MPa/yr (solid lines), using 
AS = 0.1 MPa and Ac7=0.01 MPa. Colors indicate the stressing rate. Steady-state 
seismicity rate r for each stressing rate is also shown (dashed lines). As stressing rate 
increases, R increases while fa decreases, b) Cumulative number of events over time 
obtained by integrating curves in Fig. 2a. The difference between solid and dashed lines 
of similar color represents aftershocks due to the sudden stress change that is present in 
the solid curves, c) Number of aftershocks N produced by the uniform stress change at 
each of the relative stressing rate values S/%, relative to number % produced at the 
background stressing rate. The least-squares fit shows that as stressing rate increases, 
the number of aftershocks also increases, indicating that the parameter K is dependent 
on stressing rate. 


seismicity rate R following a uniform stress step at different magnitudes 
of stressing rate relative to some background stressing rate Sb, ranging 
from 1 to 1000 (Fig. 2a, solid lines). In agreement with the first 
simulation, the results show that as the stressing rate increases, the 
seismicity rate increases, while ta decreases. 


(a) 
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2.3. Combining the ETAS and rate-state models 

To combine the ETAS and rate-state models into a single model 
appropriately, we now examine the relationships between parameters 
of both models and their dependence on stressing rate. The 
parameters we will consider are the ETAS parameters a, p, c, K, and 
p. The dependence of these parameters on stressing rate will 
determine the dependence of Rc on stressing rate, which in turn will 
determine the functional form of R that we seek to establish. 

The ETAS parameter a describes the efficiency of an earthquake of 
a given magnitude at generating aftershocks. It has no direct 
equivalent in the rate-state model, which incorporates no magnitude 
dependence in its equations. However, there is an implicit magnitude 
dependence in the rate-state model, in that larger earthquakes 
increase the stress-levels in a greater volume of the crust than small 
ones. Therefore, we assume that a is related to the spatial extent of a 
stress step and independent of stressing rate. 

Simulations using Eq. (5) show that the ETAS parameter p is 
essentially 1 and independent of stressing rate in the rate-state 
model. Both theoretical and observational studies also suggest that 
this Omori decay parameter is more influenced by factors that are 
relatively stressing rate independent, such as heterogeneity in 
temperature/heat flow or structure (e.g., Mogi, 1962, 1967; 
Kisslinger and Jones, 1991; Utsu et al., 1995). Recently, Helmstetter 
and Shaw (2006) have also shown thatp can be related to the rate- 


state parameter Acj and the spatial distribution of the stress field on 
a fault. Therefore, as p seems to be more sensitive to longer-term 
heterogeneities on a fault, in our model we consider p independent 
of stressing rate. 

In the rate-state model, the ETAS parameter c can be analytically 
related to rate-state parameters and stressing rate (Dieterich, 1994). 
In reality, it is difficult to obtain an accurate measurement of c because 
of tbe incomplete detection of early aftershocks (Utsu et al., 1995). 
Therefore, any dependence that c may have on stressing rate will most 
likely be obscured by this effect, and so we consider c independent of 
stressing rate. 

The last two ETAS parameters (K and p) do not have as 
straightforward a relationship with stressing rate. The rate-state 
model predicts that background seismicity rate (i.e., seismicity not 
triggered by an earthquake) depends on stressing rate. As the 
stressing rate increases, so does the background seismicity rate 
(blue lines in Fig. 1, dashed lines in Fig. 2). The ETAS model on the 
other hand assumes that background seismicity rate p is constant in a 
particular time interval. 

The ETAS model also assumes that aftershock productivity K is 
independent of stressing rate. On the other hand, the rate-state model 
predicts that K increases with stressing rate. Therefore, an earthquake 
with a given stress drop that occurs during a time of lower stressing 
rate will produce fewer aftershocks than if it had occurred during a 
time of higher stressing rate. 
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Fig. 3. Maps of seismicity used in analysis, a) M> 1.9 events in the Obsidian Buttes region from 1985-2005. Events in the 2005 swarm are shown in magenta, b) M>1.5 events in the 
Kilauea region from 2001-2007. Events in the 2005 swarm are shown in magenta, c) M>2 events in the Boso region from 1992 to 2007. Events in the 2002 swarm are shown in 
magenta, events in the 2007 swarm are shown in cyan. 
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The simulations using Eq. (5) demonstrate this behavior. The 
seismicity rates associated with the five different stressing rates in 
Fig. 2a are integrated to determine the cumulative number of events 
for each stressing rate (Fig. 2b). When comparing the case with a 
stress step (solid lines) to the case with no stress step (dashed lines), 
the difference between the two curves at a given stressing rate is due 
to the aftershocks produced by the stress step. We can then compare 
the number of aftershocks N produced by uniform stress changes at 
the five different stressing rates to the number of aftershocks Nb 
produced by the same stress change at the background stressing rate 
(Fig. 2c). As tbe stressing rate increases, the number of aftershocks 
produced also increases. For example, an earthquake that occurs at a 
stressing rate 1000 times higher than the background stressing rate 
produces -700 times more aftershocks than a similar sized earth¬ 
quake occurring at the background stressing rate when 
Aa= 0.01 MPa. The ratio N/Nb thus reflects the increase in aftershock 
productivity K predicted by the rate-state model due to the increase 
in stressing rate unlike the ETAS model, in which K is not related to 
stressing rate. Therefore, the main issues that need to be resolved in 
order to build a consistent combined model of seismicity rate involves 
ascertaining the dependence of K and background seismicity rate p on 
stressing rate. 

3. Data analysis 

Many studies suggest that swarms are a response to geodetically 
observed increases in stressing rate (e.g., Lohman and McGuire, 2007; 
Ozawa et al., 2007). Therefore, by analyzing swarms, we can establish 
whether the ETAS parameters I< and p change during periods of high 
stressing rates. We use three types of analyses: first, we fit the ETAS 
model to a catalog containing a swarm to determine if the triggering 
behavior is non-stationary during swarms. Second, we divide the 
catalog into pre-swarm and swarm portions, fit the ETAS model to 
each and compare the parameter estimates to determine how they 
change during swarms. Finally, we compare aftershock counts of a 
moderate-sized earthquake that occurred during a stressing rate 


transient to aftershock counts of other earthquakes in the catalog to 
test the rate-state model prediction that aftershock productivity K 
scales with stressing rate. 

We examine 4 earthquake swarms that geodetic studies have 
linked to changes in stressing rate: the 2005 Obsidian Buttes swarm 
(Lohman and McGuire, 2007), the 2005 Kilauea swarm (Wolfe et al., 
2007), and the 2002 and 2007 Boso swarms (Ozawa et al., 2003, 
2007). Catalogs for these swarms were obtained from the Southern 
California Earthquake Data Center, the Advanced National Seismic 
System, and the Japan Meteorological Agency respectively. 

3.1. Detection of anomalous seismicity rates 

The ETAS model when used as a diagnostic tool can identify time 
periods when seismicity rates do not behave as typical aftershock 
sequences (Ogata, 1988; McGuire et al., 2005). We apply the method 
described in Ogata (2005) by fitting the ETAS model to a catalog that 
includes a swarm. We then employ a transformation described in 
Ogata (1988) which utilizes the following theoretical cumulative 
function: 

A(t)= /'R(s)ds (7) 

.10 

where R is the seismicity rate predicted by the ETAS model (Eq. (1)). 
The occurrence times q in the catalog are transformed into Tj = A(ti). If 
the earthquakes in the catalog are described well by the ETAS model, 
then q will be distributed according to a stationary Poisson process, 
and a plot of the actual cumulative number of events vs. the 
theoretical number of events (i.e., the transformed time r) will be 
linear. Anomalous seismicity that the ETAS model cannot explain will 
appear as deviations from this trend. 

3.1.1. 2005 Obsidian Buttes swarm 

In August 2005, an earthquake swarm occurred over the course of 
two weeks on a continental strike-slip fault in the Salton Trough in 
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Fig. 4. Results of optimization of the ETAS model for the 2005 Obsidian Buttes catalog, comparing the cumulative number of events to the transformed time (ETAS predicted number 
of events). The observed data are shown in blue and the ETAS prediction in red. Bottom panels show the magnitudes of the events in the swarm. The ETAS model is optimized until 
just prior to the swarm and extrapolated for the remainder of the catalog. Black lines signify the 2a bounds of the extrapolation. A significant deviation from the ETAS prediction 
occurs near the beginning of the swarm, (for interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 
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Transformed Time 


Fig. 5. Results of ETAS model optimization for the 2005 Kilauea catalog. See Fig. 4 for 
symbol explanation. The ETAS model ceases to adequately fit the catalog early in the 
swarm. 


southern California. Lohman and McGuire (2007) concluded that the 
swarm was triggered primarily hy aseismic fault creep that released 
moment equivalent to a 5.7 earthquake and increased the 
stressing rate (and seismicity rate) hy about a factor of 1000. 


The earthquake catalog used in our analysis consists of events from 
1985-2005, with a magnitude of completeness M(.= 1.9 (Fig. 3a). The 
ETAS model was optimized through 2005 until just prior to the swarm 
and then extrapolated until 2006. The transformed time plot shows 
that a significant deviation from the ETAS predicted trend occurs at 
the time of the swarm (Eig. 4). More events occurred during the 
swarm than the ETAS model can explain with the parameters that best 
fit the preceding catalog. This suggests that at least one of the ETAS 
parameters changes during periods of high stressing rate. 

3.1.2. 2005 Kilauea swarm 

Slow earthquakes that trigger microseismicity periodically occur 
on the south flank of Kilauea Volcano in Hawaii (Cervelli et al., 2002; 
Brooks et al., 2006; Segall et al., 2006; Wolfe et al., 2007). In this study, 
we focus on a slow earthquake that occurred on 26 January 2005 and 
released moment equivalent to a Mw 5.8 earthquake over the course 
of hours to days (Brooks et al., 2006; Wolfe et al., 2007). 

The catalog we analyze consists of events occurring from 2001- 
2007, with a catalog completeness of Me = 1.5 (Fig. 3b). We optimize 
the ETAS model from 2001-2005 and extrapolate through the 
remainder of the catalog. Again, the results show that a significant 
deviation from typical aftershock behavior occurs at the time of the 
swarm (blue curve above the black confidence limits in Eig. 5), 
suggesting a stressing rate dependence of one or more parameters. 

3.1.3. 2002 and 2007 Boso swarms 

The Boso peninsula in central japan has been the site of recurring 
slow slip events on the subduction thrust interface in 1996, 2002 and 
2007 (Ozawa et al., 2003; Sagiya, 2004; Ozawa et al., 2007). These 
events, detected by GPS instruments, lasted on the order of a week and 
were accompanied by earthquake swarm activity. Ozawa et al. (2007) 
suggest that the slow slip events are the primary driving process 
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Fig. 6. Results of ETAS model optimization for the 2002 and 2007 Boso swarms. See Fig. 4 for symbol explanation. For both swarms, the ETAS model again fails to explain the amount 
of seismicity that occurs during the swarms, indicating that the ETAS parameters that best fit the catalog prior to the swarms no longer fit during the swarms. 
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Table 1 

Comparison of maximum likelihood estimates of ETAS parameters before and during 
each swarm. 



ETAS model is optimized from 1992 to February 2007 and extrapolated 
through 2008. Due to the short duration of the 2002 swarm, it should 
have very little effect on the parameter estimates. The results indicate 
that anomalous seismicity rates occur during the slow slip events in 2002 
and 2007 that cannot be explained by the ETAS model (Fig. 6). Again, this 
suggests that at least one ETAS parameter depends on stressing rate. 

3.2. Fitting ETAS to earthquake swarms 


triggering the swarms, similar to the Obsidian Buttes swarm (Lohman 
and McGuire, 2007). 

Our catalog consists of M>2 events from 1992 to 2007 (Fig. 3c). To 
obtain the best-fitting parameter estimates for the catalog as a whole, the 


One way to determine which ETAS parameters change during 
swarms (i.e., high stressing rate periods) is to fit the ETAS model to the 
pre-swarm portion of the catalog and compare it to ETAS fit to the 
swarm alone. Table 1 shows the pre-swarm and swarm estimates of 
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Fig. 7. Results of applying the ETAS model only on swarm seismicity from a) the 2005 Obsidian Buttes catalog, b) the 2005 Kilauea catalog, c) the 2002 Boso catalog, and d) the 2007 
Boso catalog. In all 4 cases, the ETAS model requires increases in K and compared to pre-swarm estimates in order to adequately fit the data during the swarm. 
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the ETAS parameters for each swarm. In most cases, in order for the 
model to converge, p was held fixed at 1.0. The observed and predicted 
cumulative numbers of events for each swarm are shown in Fig. 7a-d. 
The poorer quality of the fits could suggest that a time-dependent p 
may be necessary to more accurately fit the data. For all of the swarms, 
the ETAS model finds changes in K by factors of 2-4. Flowever, the 
parameter p increases by 1-3 orders of magnitude during the swarms. 
Therefore, with the ETAS model, stressing rate transients appear to 
primarily increase the background seismicity rate without increasing 
aftershock productivity substantially. 

3.3. Comparison of rate-state predictions with observations 

A final way to test the stressing rate dependence of K is to look at 
moderate sized earthquakes that occur during a swarm. By counting 
the number of aftershocks following these earthquakes (in narrow 
space-time windows) and comparing to the number of aftershocks 
produced by other earthquakes in the catalog (during low stressing 
rate periods), we can test the hypothesis that K is stressing rate 
dependent. Assuming that the ETAS parameters a, p, and c remain 
constant over time, the average number of aftershocks following an 
earthquake of magnitude M can be expressed as: 

N= -^10“<“-“'> (8) 

1 — n ' ' 

where the branching ratio n=Kb/(b — cx) (Helmstetter and Sornette, 
2003; McGuire et al., 2005). Therefore, if the ETAS and Gutenberg- 
Richter parameters are assumed to remain constant throughout the 
catalog, N is primarily a function of the difference between mainshock 
magnitude M and catalog completeness threshold Me, and a plot of the 
logarithm of aftershock counts of mainshocks in the catalog will be linear 
with respect to the mainshock magnitudes (McGuire et al., 2005). In 
contrast, the rate-state equations predict a greater productivity (larger 
K) during the transient and Eq. (8) will not describe the data well. 

Fig. 8 shows the aftershock productivity for the Obsidian Buttes 
catalog. Aftershocks in a 1-day time window were counted for 



Fig. 8. Aftershocks per mainshock vs. difference between mainshock magnitude and 
magnitude of completeness Me =1.9 for events of Mmain^4 in the Obsidian Buttes 
catalog from 1985-2005 (circles). Aftershocks for each event were counted within a 
1 day window. Lines depict lines of constant aftershock parameters, where number of 
aftershocks depends on Mmain~^c- Solid line shows least-squares fit to data; dashed 
line shows increase in aftershock productivity predicted by rate-state model for a 
stressing rate of 1000 times background stressing rate. The actual number of 
aftershocks following the M5.1 mainshock (star) is much less than the number 
predicted by the rate-state model (square) and falls on a line consistent with other 
mainshocks in the catalog, suggesting that K (i.e., aftershock productivity) is not 
stressing rate dependent. 


mainshocks with M>4, occurring sufficiently apart in time so as not 
to interact with one another. During the swarm, Lohman and McGuire 
(2007) estimated that a stressing rate transient of -1000 times the 
background stressing rate occurred. Therefore, the rate-state model 
equations predict that the M5.1 earthquake that occurred during the 
swarm should produce almost 1000 times more aftershocks than a 
similar sized earthquake occurring at typical stressing rates (Fig. 2c). 
Flowever, the actual number of aftershocks observed following the 
earthquake was not that large (star in Fig. 8). The aftershock count for 
this event in fact plots on the same constant line as the other events in 
the catalog, suggesting that K is independent of stressing rate. A 
concern is that the lack of increase in K could be due to the incomplete 
detection of early aftershocks. However, we have carefully taken the 
magnitude of completeness Me into account for each of the catalogs in 
our analysis. Moreover, the rate-state model equations and thus 
predictions are defined for a given magnitude interval that we assume 
to be M>Me (Dieterich, 1994). Therefore, undetected aftershocks are 
unlikely to be the primary reason for the lack of an increase in K. 

3.4. Summary 

We have analyzed 4 different earthquake swarms to examine the 
dependence of the ETAS parameters K and p on stressing rate. The 
ETAS model identified the swarms as anomalous seismicity that 
cannot be fit with the same parameters as the rest of the catalog, 
suggesting that at least one of the parameters changes with stressing 
rate. However, when the ETAS model was fit to the swarms alone, 
estimates for K changed very little compared to the pre-swarm fit 
while the estimates for p increased by several orders of magnitude. 
Finally, the aftershock count following the M5.1 Obsidian Buttes 
earthquake revealed no substantial increase in K during the 
heightened stressing rate associated with the swarm. Together these 
results suggest that stressing rate transients increase the background 
seismicity rate p without causing a substantial increase in aftershock 
productivity K. 

4. Discussion and conclusion 

The primary conflict between the ETAS and rate-state models of 
seismicity rate lies in the dependence of aftershock productivity on 
stressing rate. Our results suggest that, contrary to rate-state model 
predictions, the aftershock productivity is unaffected during stressing 
rate transients, which instead increase the background seismicity rate. 
The key to this discrepancy lies in the evolution of the rate-state 
variable y. In the rate-state model, an increase in seismicity rate 
implies that y has evolved to a new steady-state value (Eq. (2)). Our 
simulations show that with this increase in seismicity rate comes an 
increase in aftershock productivity (Eig. 2). Therefore, the lack of an 
increase in productivity suggests that the state variable y has not 
evolved. Thus, there is a fundamental conflict between the heightened 
seismicity rate and the unchanged aftershock productivity we have 
observed in our analysis of swarms. 

An additional complication is that the evolution of y also depends 
on the frictional parameter Act (see Eq. (3)). This parameter controls 
how quickly y evolves in response to changes in stressing rate, and 
therefore ultimately affects the aftershock productivity of a stress step. 
We can again use Eqs. (2)-(4) to explore in detail how the change in 
aftershock productivity with stressing rate varies with Aa. We 
compare two stress histories, one in which an earthquake 
(AS = 1 MPa) occurs during a background stressing rate of 0.2 MPa/ 
yr, and one in which a similar stress step occurs three days after a 
stressing rate transient begins. We use Eqs. (3) and (4) to calculate y 
for each stress history and Eq. (2) to obtain seismicity rates that can 
then be integrated to estimate the number of aftershocks produced by 
each earthquake. Fig. 9 compares the number of aftershocks N 2 
produced by an earthquake that occurs three days after a stressing rate 
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Fig. 9. Change in aftershock productivity N 2 /N 1 vs. Aa following various magnitudes of 
stressing rate jumps relative to the background stressing rate of 0.2 MPa/yr (symbols), 
for an earthquake with AS = 1 MPa. The ratio predicted by the ETAS fits in Table 1 for the 
Obsidian Buttes swarm is indicated. Values for Aa are also shown, given laboratory 
values of A for quartz (Chester and Higgs, 1992) and granite (Blanpied et al., 1998) in 
hydrostatic conditions at a depth of 4 km for temperatures ranging from 300 °C to 
600 °C. For any given range of Aa, the rate-state model cannot satisfy both the increase 
in seismicity/stressing rate and the lack of change in aftershock productivity observed 
during the Obsidian Buttes swarm. 


transient begins to the number of aftershocks Nj produced by an 
earthquake that occurs during the background stressing rate, using 
different values of Aa ranging from 10^^ to 3 MPa. The ratio N 2 /N 1 
essentially gives the expected increase in aftershock productivity K 
during a stressing rate transient. Fig. 9 demonstrates that the 
predicted change in aftershock productivity is highly dependent on 
the value of Aa used in the rate-state equations. 

Catalli et al. (2008) recently examined the role of Atrin modeling 
seismicity rate variations and found that it controlled the total number 
of aftershocks triggered by an earthquake primarily in two ways. First, 
Acj controls the instantaneous change in 7 (and therefore seismicity 
rate) due to a sudden stress step (Eq. (4)), so that as Acj increases, the 
instantaneous change in seismicity rate decreases. Second, the 
duration of aftershock sequences, ta, also depends on Act; as Aa 
increases, ta increases (Eq. ( 6 )). The simulations in this study demon¬ 
strate that these two effects are also dependent on stressing rate 
(figs. 1-2). As stressing rate increases, the change in seismicity rate 
due to a stress step increases, but the aftershock duration ta decreases. 
Therefore, the range of aftershock productivity behavior seen in fig. 9 
reflects the tradeoffs in how these two effects are controlled by both 
Aa and stressing rate. 

To compare this predicted behavior with a real-life example, for 
the 2005 Obsidian Buttes M5.1 earthquake, which occurred three days 
after an increase in stressing rate of almost three orders of magnitude, 
we found N 2 /N]~l from aftershock counts (fig. 8 ). Additionally, the 
ETAS model fitting resulted in an increase in background seismicity 
rate by over three orders of magnitude (Table 1), which agrees with 
the observed increase in seismicity rate. Fig. 9 shows that the rate- 
state model cannot satisfy all of these observations in a small range of 
Acj. Typical estimates of Acj from earthquake catalogs range from 10^^ 
to 10^’ MPa (e. g.. Gross and Kisslinger, 1997; Harris and Simpson, 
1998; Toda et al., 1998; Belardinelli et al., 1999; Console et al., 2007). In 
this range of Acj, 7 evolves quickly, so that jumps in stressing rate 
cause Jumps in seismicity rate, but also cause jumps in N 2 /N 1 (i.e., 
aftershock productivity). Laboratory measurements of A for quartz 
and granite (Chester and Higgs, 1992; Blanpied et al., 1998) result in 
higher values of Acj (10^ ’ to 1 MPa) for faults under hydrostatic pore 
pressure at a depth of 4 km. At these values of Acj, although aftershock 
productivity does not change with the jump in stressing rate, neither 


does the seismicity rate, because 7 has not evolved to any great extent. 
We find that the increase in seismicity rate and the lack of change in 
aftershock productivity observed in the Obsidian Buttes swarm cannot 
both be satisfied simultaneously using the rate-state model, because 
the two observations imply fundamentally different things about 
whether 7 has evolved or not. Therefore, some caution is necessary 
when applying the rate-state inversion algorithm to obtain stressing- 
rate changes from earthquake catalogs. 

Given our observations of the dependence of the ETAS parameters 
on stressing rate, we can now specify a combined ETAS/rate-state 
model of seismicity rate to detect stressing rate transients from 
earthquake catalogs. As described earlier, the seismicity rate R in a 
catalog is a function of an aseismically-triggered component and an 
earthquake-earthquake triggered component Rc. While Ra is clearly 
related to stressing rate, the relationship between Rc and stressing 
rate was unclear. Our results suggest that K is independent of stressing 
rate for a particular region, and so Rc is independent of stressing rate. 
R can then essentially be separated into the aseismic component Ra 
and the coseismic component Rc (represented by the ETAS model); 

R = R.+R. = R.+ y" 7 -CFT (9) 

Then Ra is effectively a time dependent version of the ETAS 
parameter p (see Eq. (1)), and to obtain it, one can simply subtract the 
ETAS-predicted Rc from the observed rate R. The residual Ra can then 
be directly related to a stressing rate Sa caused by aseismic 
deformation through the rate-state model equations; 


Ra — R — Rc — k — 
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(t-ti + cf 


r 

5r7 


( 10 ) 


dy = 


dt 

Aa 


[1 - t(Sa 



( 11 ) 


where Sb is the background tectonic stressing rate. The use of the ETAS 
model to estimate Rc reduces the impact of aftershock sequences on 
the estimation of the aseismically-triggered seismicity rate, while the 
rate-state model establishes the relationship between aseismically- 
triggered seismicity rates and stressing rates. 

There are a number of caveats to keep in mind about this model. 
First, a fundamental assumption this model makes is that the ETAS 
parameters R, a, c and p are constant in space and time and can 
describe all coseismically-triggered seismicity in a catalog completely. 
These parameters can in fact vary from sequence to sequence, as well 
as place to place (Ogata, 1998). However, our data analysis suggests 
that at least within relatively small and homogeneous regions, these 
parameters will remain constant, and therefore changes in the 
observed seismicity rate will be primarily mapped into changes in 
Ra. Second, practical applications of this model will have to be careful 
about catalog completeness, to ensure that undetected events such as 
early aftershocks do not affect the results. An additional point to bear 
in mind is the need to smooth the seismicity rates over some time 
window. Currently there is no way to estimate Ra at arbitrarily fine 
time scales, so algorithms will need to be developed to smooth these 
estimates in an optimal way. These issues will need to be considered 
when this model is implemented in an algorithm to invert for 
stressing rate variations, but the studies of Toda et al. (2002) and 
Lohman and McGuire (2007) suggest that our model is correct at least 
at the order of magnitude level. To verify the rate-state equations for 
background rate at a higher level of precision would require high 
sample-rate geodetic measurements (e.g., strain- or tiltmeters). 
However, to first order, our model provides a simple and direct way 
to quantitatively relate aseismic stressing rate transients to seismicity 
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data and will allow future studies to invert seismicity catalogs to 
detect stressing rate variations caused by transient aseismic processes. 
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Chapter 4: Detecting Aseismic Strain Transients From Seismicity Data 


Abstract 

Aseismic deformation transients such as fluid flow, magma migration, and slow 
slip can trigger changes in seismicity rate [e.g., Dieterich et al, 2000; Toda et al, 2002; 
Hainzl and Ogata, 2005; Lohman and McGuire, 2007; Ozawa et al., 2007], W e present a 
method that can detect these seismicity rate variations and associate these anomalies with 
underlying variations in stressing rate. Because ordinary aftershock sequences often 
obscure changes in the background seismicity caused by aseismic processes, we combine 
the stochastic Epidemic Type Aftershock Sequence model [Ogata, 1988; Ogata, 1998] 
that describes aftershock sequences well and the physically based rate- and state- 
dependent friction seismicity model [Dieterich, 1994] into a single seismicity rate model 
that models both aftershock activity and changes in background seismicity rate. We 
implement this model into a data assimilation algorithm that inverts seismicity catalogs to 
estimate space-time variations in stressing rate. We evaluate the method using a 
synthetic catalog, and then apply it to a catalog of M>1.5 events that occur in the Salton 
Trough from 1990-2009. We validate our stressing rate estimates by comparing them to 
estimates from a geodetically-derived slip model for a large creep event on the Obsidian 
Buttes fault [Lohman and McGuire, 2007]. The results demonstrate that our approach 
can identify large aseismic deformation transients in a multi-decade long earthquake 
catalog and roughly constrain the absolute magnitude of the stressing rate transients. 
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1. Introduction 

Transient aseismic processes such as fluid flow, magmatic intrusions, or slow slip 
events can alter the stress state and trigger seismicity in a variety of tectonic 
environments. For example, earthquake swarms associated with aseismic deformation 
have been observed in subduction zones such as Japan and New Zealand [Ozawa et al, 
2007; Delahaye et al., 2009], continental strike-slip faults [Lohman and McGuire, 2007], 
and volcanic regions [Dieterich et al., 2000; Toda et al., 2002; Segall et al., 2006; Wolfe 
et al., 2007; Montgomery-Brown et al, 2009]. However, there is not a straightforward 
relationship between the magnitude of the aseismic transient and the amount of seismicity 
triggered. In the Salton Trough of California, a M w 5.7 creep event triggered a swarm of 
-1000 M <5.1 earthquakes [Lohman and McGuire, 2007], whereas offshore of central 
Honshu, japan, recurrent M w~6.5 slow earthquakes on the plate interface typically trigger 
swarms of -lOs of M <5 earthquakes [Ozawa et al, 2007]. Moreover, a Mw 4.7 creep 
event detected on the Superstition Hills fault in southern California did not trigger any 
earthquakes [tye/era/., 2009]. There is evidence, however, that the seismicity rate varies 
in proportion to the stressing rate increase caused by the aseismic transient [Toda et al, 
2002; Segall et al, 2006; Lohman and McGuire, 2007]. Although this is likely 
complicated by a dependence on ambient stress conditions such as the timing within a 
fault's seismic cycle, there appears to be a possibility of utilizing seismicity rate 
variations observed in earthquake catalogs to detect when and where aseismic processes 
are occurring on some faults, which would ordinarily require high-quality geodetic data 
to resolve. Moreover, the variations in earthquake triggering observed during aseismic 
deformation transients suggest an opportunity to learn about the physics of earthquake 
triggering, and to constrain the extent to which seismicity is triggered by potentially 
observable geophysical processes, such as slow slip, and hence may be forecastable in a 
probabilistic sense. 

The rate- and state-dependent friction model [Dieterich, 1994; Dieterich et al, 
2000] is one approach for mapping seismicity rate variations to underlying stressing rate 
variations. This model has been used to estimate stress changes caused by a dike 
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intrusion [Dieterich et ah, 2000] and slow slip events [Segall et ah, 2006; Montgomery- 
Brown et al, 2009] on Kilauea volcano as well as a large intrusion in the Izu islands 
[Toda et al, 2002], Oftentimes, however, the background seismicity rate changes due to 
aseismic processes can be obscured by the aftershock sequences in a catalog, particularly 
in regions such as southern California that are characterized by high aftershock 
productivity [Helmstetter, 2003; Helmstetter and Sornette, 2003a]. In SUCh cases, the 
stressing rate changes estimated by the rate-state inversion will be a combination of those 
due to the underlying aseismic process and those due to the static stress changes resulting 
from the seismicity itself. Moreover, the rate-state equations predict that changes in 
stressing rate due to aseismic processes should cause many orders of magnitude changes 
in aftershock productivity, but this was not observed in earthquake swarms triggered by 
aseismic fault slip in a variety of tectonic environments [Llenos et al, 2009]. For these 
reasons, the rate-state inversion should be used with caution. 

The stochastic EpidemicTypeAftershock Sequence(ETAS) model [Ogata, 1988] 
is an alternative approach used to estimate the time-dependence of underlying driving 
mechanisms. Earthquake occurrence is modeled as a point process and characterized by 
a background rate and a set of parameters associated with Omori's law that can be 
optimized to fit a particular catalog. The ETAS model is effective at detecting when 
anomalous seismicity is being triggered by an external process [Hainzl and Ogata, 2005; 
Ogata, 2004; Ogata, 2005; McGuire etal, 2005]. However, ETAS lacks a procedure for 
estimating smooth variations in the background earthquake rate and a quantitative way to 
relate these directly to stressing rate variations. 

The rate-state model and the ETAS model have complementary strengths which, 
when combined, can provide an effective tool for detecting seismicity rate anomalies and 
relating them to underlying stress fields. In a previous study, we combined these two 
models into a single seismicity rate model that explains both aftershock activity as well as 
variations in background seismicity occurring in a region [Llenos et al, 2009]. In this 
model, the observed seismicity ratei? in a catalog is approximately a linear combination 
of an aseismically-triggered component (reflecting seismicity triggered by long-term 
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tectonic loading and aseismic processes) and an earthquake-earthquake triggered 
component (reflecting aftershock sequences). This is in contrast to the pure rate-state 
model where aftershock productivity is proportional to the stressing rate, which was not 
observed in our previous study of slow slip events, preventing us from using it for the 
coseismically-triggered component [Llenos et al, 2009]. We instead use the ETAS 
model to estimate the coseismically-triggered component of seismicity rate and subtract it 
from the observed (binned) catalog rate. We then relate the residual seismicity rate 
directly to an aseismic stressing rate using the rate-state model. 

We implement this combined model in a data assimilation algorithm that uses 
seismicity catalogs to estimate stressing rate variations. We incorporate the rate-state 
equations into a state-space model and use an extended Kalman filter to estimate spatial 
and temporal variations in stressing rate. We evaluate this method with a synthetic test, 
and apply it to a catalog of events from the Salton Trough in southern California in which 
a geodetically-detected aseismic transient occurred. To validate our approach, we 
compare the estimated stressing rate peak during this aseismic slip event to that obtained 
from a slip inversion of geodetic data [Lohman and McGuire, 2007]. The results suggest 
that our seismicity inversion method provides an accurate way to detect and locate 
transient deformation strictly from seismicity catalogs and may provide a means to 
constrain the absolute magnitude of stressing rate variations. 

2. Method 

The seismicity rate 7? observed in a catalog is in general a function of the stressing 
rate acting on the population of faults in a region [Dieterich, 1994]. We assume that the 
seismicity in a catalog is primarily triggered by three mechanisms: earthquake-earthquake 
interactions (i.e., as aftershocks), transient changes to the stress field from aseismic 
processes such as fluid flow or slow slip, and long-term tectonic loading. Therefore, 
R{t,x) =f{S) =f{R^,R,^,Rj.), where Ra reflects seismicity triggered by aseismic 

processes, Rc represents a coseismically-triggered component (i.e., earthquakes triggered 
by other earthquakes), and is triggered by long-term tectonic loading. We assume that 
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Rt is small relative to Ra if an aseismic deformation/fluid-flow transient is occurring and 
therefore combine it with Ra, so that R{t,x) ~f{R^,R^). In previous studies, the ETAS 
model and the rate-state model have each been used to investigate both Ra and Rc- 


2.1 ETAS model 

The space-time ETAS model is a point process model that represents the 
seismicity rate R{t,x,y) as a summation of the aftershock sequences produced by each 
event prior to time t plus a time-independent background seismicity rate ^[Ogata, 1998; 
Zhuang et al, 2005; Ogata and Zhuang, 2006]: 

R [t, X, y ) =/j{x, y ) + ^ k[m . -t^)- f [x -x. .y-y^M^) ,, > 



( 2 ) 

(3) 

(4) 


The function a(M) is the expected number of events triggered by an earthquake of 
magnitude M and involves the aftershock productivity K and scaling parameter a, the 
function ifj{t) is the probability density function form of the modified Omori law and is 
specified by the Omori decay parameters c and p [Omori, 1894; Utsu, 1961], and the 
function f{x,y;M) describes the spatial distribution of the events triggered by an 
earthquake of magnitude M, specified by the parameters d, q, and scaling parameter q. 
The ETAS parameters K, c, p, a, d, q, and qate generally obtained using maximum 
likelihood estimation from the observed occurrence times u and magnitudes M, of 
earthquakes in a catalog with a magnitude of completeness Me [Ogata and Zhuang, 
2006], 
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2.2 Rate- and state-dependent friction model 

In the rate- and state-dependent friction model, the seismicity ratei? is a function 
of the stressing rate s acting on a population of faults governed by rate- and state- 
dependent friction [Dieterich, 1994]: 



dr=4^k-s] ( 6 ) 

Ao 

where r is a steady-state reference seismicity rate associated with a reference stressing 
rate 5^, 5 is the modified Coulomb stressing rate, yls a state variable, and A is a fault 
constitutive parameter. We assume the normal stress a and therefore the frictional 
parameter Aaremain constant. 


2.3 Combined ETA S/rate-state model 

By examining how the ETAS and rate-state model parameters changed during 
periods of high stressing rate, Llenos et al. [2009] determined that aseismic transients 
primarily increase the background seismicity rate (i.e., ETAS parameter /J) while other 
aftershock parameters such as productivity K remain relatively unaffected by the increase 
in stressing rate. The seismicity ratei? can then be (approximately) modeled as a linear 
combination of the aseismically-triggered component Ra and the coseismically-triggered 
component Rc‘. 


R{t)=R,{t)+R,{t)=^R,{t)+Y, 

tj S 


{t -t. +cy 


(7) 


where Rc is represented by the ETAS model, and Ra is essentially a time-dependent 
version of /ithat will reflect the variations in seismicity rate triggered by transients. We 
approximate Ra by subtracting out an ETAS-predicted rate from the observed rate i? in a 
catalog. We then directly relate our estimate of Ra to an aseismic stressing rate 
through the rate-state model equations. We utilize the space-time version of the ETAS 
model [Ogata andZhuang, 2006], SO that the model becomes: 
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( 8 ) 

c J rfe"'"''‘■'^' j s^Y 

dY=^V-\{s,+sJ (9) 

Aa 

where is the background plate tectonic stressing rate. This combined model of 

seismicity rate now reduces the impact of aftershock sequences on the estimation of Ra 
by using ETAS to model Rc and provides a way to quantitatively relate seismicity rate 
variations to stressing rate variations. 


2.4 Extended Kalman filter algorithm 

We incorporate the aseismically-triggered ratei?A into a state-space model, which 
describes a system using a state vector that evolves over time and can be related to 
observable data. Our state-space model consists of a non-linear measurement equation 
that relates the observations d*: to the state vector x* at each ti me step k\ 

^k=KM+H’ (4 ~n{0,rJ ( 10 ) 

and a non-linear state transition equation that describes how the state vector evolves over 
time: 

(**: ) ^k-A ~ (H) 

Our State variables Xk consist of the long-term tectonic loading rate s^, the aseismic 

stressing rate , and the logarithm of the rate-state seismicity state variable )/to ensure 
positivity. We also divide the region up into N spatial boxes and estimate the state 
variables in each one, so that the full state vector is: 

X. =1*^,, k I ih linn k I k I ih )- 1 n K2 ih I ■ ■ 

s,^tkls^^hh^VNih)] (12) 

The variables and ln( j) are modeled as random walk processes with scale parameters 
of rand £ respectively. The data vector d*: consists of the aseismically-triggered ratei?A 
estimated in each spatial box and is related to the state variables through the measurement 
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matrix h/t, which is based on the rate-state model (Eq. 5). The measurement and process 
noises ujk and dk+i are assumed to be Gaussian and described by the data and process 
covariance matrices R^: and Qk+i respectively. The data covariance matrix Rt is assumed 
to be known only up to a scalar multiplier such that , where It is a matrix 

that contains any known correlations between the spatial cells and components [Segall 
and Matthews, 1997], I n our case, we assume ~1 and estimate Ik a priori based on the 
variance of the data vector during stable time periods. 

The system is solved using an extended Kalman filter algorithm [e.g., Gelb, 1974; 
Anderson and Moore, 2005]. We first make a priori estimates of the State vector X 0|0 and 
covariance Io|o and update them with the first observations to obtain xi|i and Ii|i. We 
then use the filter prediction equations to estimate the state vector and covariance at the 
next time step k+l\ 

2.^1. (14) 


where T is the Jacobian of the non-linear state transition model t, which incorporates the 
rate-state evolution equation (Eq. 6): 

T, (15) 

dx ^ 

and the elements of the 3W x3W process covariance matrix Q/t+i are determined by the 
properties of the random walks used to model the state variables 5^ and ln(}): 




0 0 
0 

0 0 

0 


0 

0 


0 



(16) 


The predictions are then compared with the observations at the current time step and 
updated using the filter update equations: 

'^k-a\k-A ^U/t)i (1'^) 
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where the Kalman gain g is computed by: 


and H is the Jacobian of the non-linear measurement matrix h, which incorporates the 
rate-state model equation (Eq. 5): 


H 


k 



( 20 ) 


These prediction and update steps are iterated forward through the entire data set, 


resulting in estimates of the state variables xm at each time step k given the data up to 
time step k. The filter can be run backwards in time in a process known as smoothing to 
obtain a backsmoothed estimate of the state variables (i.e., Xkw, or estimates of x at time 
step given the entire data set). 


2.5 Likelihood calculations 

There are a total of 11 time-independent parameters in our model, seven of which 
are the space-time ETAS parameters {K, c, p, a, d, q, and /;), two of which are rate-state 
parameters (Aaand r), and two of which are filter hyper-parameters (rand £). We first 
optimize the model over the set of ETAS parameters ^given a history of occurrence Ht 
using a point process likelihood function [Ogata, 1998; Daley and Vere-Jones, 2002]: 


\nL{(p\H,)=^\n R[t,.x,. y, I H,) y | H, )dAdt (21) 

0 A 

While this likelihood computation works well to estimate the ETAS parameters, it is 
relatively insensitive to changes in the non-ETAS parameters. Based on synthetic tests, it 
is preferable to optimize the set of non-ETAS parameters ^by computing a likelihood 
based on prediction-error decomposition [Harvey, 1989; Segall and Matthews, 1997]. 
For Gaussian data, the likelihood can be expressed as: 


logL(^) -N, logArJ-^f;iog|V,| -\n, log 


t 




where ATrf is the total number of observations, 14 is the innovation at each time step k. 


( 22 ) 
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( 24 ) 


and Yk is the variance of the innovation: 
V, +R, 


This iikeiihood function aiso empioys a maximum iikeiihood estimate for the data 
covariance scaiar muitipiier 




(25) 


With these equations, we then grid-search over muitipie non-ETAS parameter sets to 
determine their maximum iikeiihood estimate, with initiai grid sizes of a factor of 10 for 
each parameter. A fundamentai assumption we have made in this initiai impiementation 
is that the data and errors are Gaussian distributed, which is likeiy incorrect in that 
earthquake statistics typicaiiy invoive non-Gaussian distributions. Future 
impiementations couid incorporate non-Gaussian error distributions by utiiizing particie 
fiiter methods [e.g., Fukuda and Johnson, 2008; Werner, 2008; Werner et ah, 2009]. 


2.6 Summary 

Our compiete aigorithm can be summarized as foilows: 

0. Determine magnitude of completeness for the catalog. 

1. Optimize the ETAS parameters for the catalog using a space-time ETAS algorithm 
[Zhuang et al, 2005; Ogata and Zhuang, 2006]. 

2. Calculate Ra to be used as the data vector for the filter by subtracting the ETAS- 
predicted rate from the observed (binned) seismicity rate in each of spatial boxes. 

3. Given a set of non-ETAS parameters, run the extended Kalman filter to obtain the time 
history of the state variables. 

4. Compute the Gaussian likelihood associated with the set of non-ETAS parameters. 

5. Repeat steps 3 and 4 as a grid-search for different parameter sets to obtain maximum 
likelihood estimates of the non-ETAS parameters. 



3. Synthetic test 

To test the validity of our algorithm, we generate a synthetic catalog of events 
using known model parameters and a stressing rate history that we will attempt to 
recover. We subdivide a 30 km by 30 km region into 9 spatial boxes (Fig. la) and 
impose a factor of 10 increase in stressing rate above a background rate of 0.1 M Pa/yr in 
the center box (Box 5). The stressing rate in the other boxes remains constant at the 
background rate (Fig. 2a). Given these stressing rate histories, we use the rate-state 
equations (Eg. 5-6) to calculate associated background seismicity rate histories, which are 
used in an ETAS simulator [Felzer et al, 2002; Helmstetter and Sornette, 2002, 2003b] 
to generate a synthetic catalog with a magnitude of completeness Me = 0. Fig. 1 shows 
the catalog event locations, magnitudes and times of occurrence. Occurrence times were 
binned in time windows of 1 day to obtain seismicity rates. The seismicity rate variations 
in both space and time due to the transient are not easily seen because the catalog is 
dominated by the aftershocks of the largest events (Fig. 2b). The anomalous seismicity 
triggered by the transient only becomes readily apparent following the removal of the 
ETAS-predicted rate (Fig. 2c). 

3.1 Parameter estimation 

We estimate the best-fitting ETAS parameters using the space-time algorithm of 
Ogata and Zhuang [2006]. The parameters used in the simulator were K = 4.5e-4 
events/day/deg^ a= 0.8, p = 1.2, c = 0.0001 days, d = 0.001 deg^ <? = 3, and /; = 0.3. 
The estimates of these parameters from the Ogata and Zhuang [2006] algorithm were.fi: = 
0.15 events/day/deg^ a= 1.8, p = 1.2, c = 0.0001 days, d = 0.5e-8 deg^ q = 1.35, and q = 
0.58. A different form for the cluster size factor tiM) (Eg. 2) is used to generate the 
synthetic catalog [Felzer et al., 2002; Helmstetter and Sornette, 2003b] than the form 
utilized in the Ogata and Zhuang [2006] algorithm, resulting in differences between the 
true and estimated values of the parameters K and a. Differences in the estimates of the 
spatial scaling parameters d, q, and q are due to the fact that the spatial probability 
density function (Eg. 4) used by the simulator also has a slightly different form than that 
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used by the Ogata and Zhuang [2006] estimation aigorithm. However, both the 
simuiation and estimation aigorithms empioyed the same temporai function (Eg. 3) based 
on the modified Omori iaw, and the true vaiues for p and c are recovered by the 
estimation aigorithm even though the assumption of constant background seismicity rate 
is vioiated. Thus, the space-time ETAS parameters appear to be weii resoived by the 
estimation aigorithm. 

The ETAS-predicted seismicity rate is subtracted from the binned seismicity rate 
to obtain the data vector (Eq. 8, 10), which is the input for the extended Kaiman fiiter. 
We ran the fiiter with a number of different vaiues for the non-ETAS parameters (r, Ao, 
T, and £) and caicuiated the Gaussian iikeiihood for each fiiter run (Eq. 22). Again, the 
fundamentai assumption made here is that the data and errors are Gaussian distributed 
(see Section 2.5). Whiie the reference seismicity rate r was weii-resoived by this grid- 
search procedure, we observed a trade-off in iikeiihood between the parameters A aand r, 
suggesting there is iittie constraint in Aa Therefore, we fixed Aato be the true vaiue 
(0.001 M Pa) and grid-searched over rand £, the scaiing parameters for the random waik 
processes used to modei and in(}) respectiveiy. 

There are two furtherconstraintsthatweempioy to seiect appropriate vaiues for r 
and £ The first constraint comes from enforcing consistency between the fiiter's 
estimates of the state variabies and in(0. Because they are modeied as stochastic 
processes with noise specified by covariance matrix Q*: (Eq. 11), they do not have to 
strictiy obey Equations 5-6. Thus, the fiiter estimate of totai stressing rate s (the sum of 
the estimates of and ) can be integrated using the rate-state equations (Eq. 5-6) to 

produce an estimate of ythat is not identicai to the fiiter estimate for y These two 
estimates of j/ideaiiy shouid cioseiy agree with one another. This piaces a constraint on 
the vaiue of r(i.e., the temporai smoothing of 5^), since it must be large enough to allow 

enough variation in S for it to integrate to y However, if ris too large (i.e., too 
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rough), the error covariance nnatrix I becomes non-positive semi-definite, and the 
prediction step (Eq. 14) can no ionger be computed. 

A second constraint is empioyed to heip choose a vaiue for £ Fig. 3 shows the 
Gaussian iikeiihood and maximum iikeiihood estimate of the data covariance muitipiier 

for various vaiues of £, caicuiated using Eq. 22 and 25. The maximum iikeiihood 

occurs when £=0.79 (Fig. 3a); however, at this vaiue of £, is significantiy iessthan 1 

(Fig. 3b). Because we estimate Rk a priori based on the variance of the data vector, 
shouid be ~ 1. Vaiues other than 1 wiii over-fit or under-fit the data. Therefore, we 

choose £such that ~1. M oreover, at higher vaiues of £ the fiiter estimate for S faiis 
to integrate to match the fiiter estimate for yas weii as at smaiier vaiues of £because the 
stressing rate estimate cannot vary enough to match the variance in K(Fig. 4a). 
Additionaiiy, higher vaiues of £uitimateiy over-fit the data (Fig. 4b). W hiie the fiiter fits 
the smaiier peaks in the data vector weii, it over-predicts the iargest peak. Smaiier vaiues 
of £ smooth out the smaiier peaks in seismicity rate but end up fitting the iargest peak 

weii. Therefore, our preferred vaiue for £is 0.16, for which <f^~l, rather than the M LE 
of 0.79. The stressing rate estimate resuits show that the estimate for s with £= 0.16 
actuaiiy comes cioser to recovering the true S than the M LE vaiue of £(Fig. 4c). 

To summarize, of the 11 modei parameters used in our aigorithm, the seven 
space-time ETAS parameters are weii-resoived by maximizing the point process 
iikeiihood function used in the Ogata and Zhuang [2006] aigorithm. The two random 
waik scaiing parameters rand £can be optimized using a Gaussian iikeiihood function 
subject to constraints. The rate-state parameters r and Aa however are not weii 
constrained by either iikeiihood function and thus shouid be fixed a priori. 

3.2 Synthetic test resuits 

Fig. 5 shows the resuiting fiiter estimates of S for each box using £= 0.16. The 
transient is cieariy identifiabie in Box 5, in both the forward and backsmoothed fiiter 


73 



estimates. Although the filter-estimated peak stressing rate is not quite as large as the 
true stress!ng rate, it is still within a factor of 2-3. The forward filter estimate of the onset 
of the transient is also delayed compared to the true onset, although the backsmoothed 
estimate mitigates this effect somewhat. Clearly, the filter is able to detect when and in 
which region the transient is occurring (Fig. 6). This is an encouraging result because the 
jump in stressing rate was only a factor of 10 while the stressing rate jumps over the plate 
tectonic rate that occur during slow earthquakes and magma intrusions are likely to be 
many orders of magnitude [Thatcher, 2001; Toda et al., 2002; Lohman and McGuire, 
2007]. 


4. Data Analysis: Salton Trough 

We next apply our algorithm to real earthquake data from the Salton Trough in 
southern California. In the Salton Trough, a transition occurs from a divergent plate 
boundary setting in the Gulf of California to the south, to the San Andreas strike-slip fault 
system to the north. The region is characterized by high heat flow [KissUnger and Jones, 
1991], which potentially acts to subdue aftershock activity [Ben-Zion and Lyakhovsky, 
2006]. A high rate of earthquake swarm activity has been observed [e.g., Richter, 1958; 
Brune and Allen, 1967; Hill et al., 1975; Johnson and Hadley, 1976; Lohman and 
McGuire, 2007; Roland and McGuire, 2009], possibly driven by magmatic intrusion 
[Hill, 1977] or aseismic fault creep [Lohman and McGuire, 2007; Roland and McGuire, 
2009]. A number of aseismic transients have also been geodetically detected in this 
region, including afterslip following the 1987 M6.6 Superstition Hills earthquake 
[Williams and Magistral , 1989], creep events on the Superstition H ills fault [Wei et al, 
2009], and aseismic creep on the Obsidian Buttes fault [Lohman and McGuire, 2007]. 

4.1 Data binninq 

We analyze a catalog of M > 1.5 earthquakes that occurred in the Salton Trough 
from February 1990 to August 2009. Because of the partial derivatives involved in 
setting up the state-space model equations (e.g., Eq. 15 and 20), we need to be able to 
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resolve the background seisnnicity rate in each space-time window. Therefore we divide 
the region up into 4 spatial boxes (Fig. 7) and bin the occurrence times into time windows 
of 20 days to obtain seismicity rates in each box. This time window is necessary to 
obtain enough earthquakes in each bin to resolve the background rate. 

4.2. Parameter estimation 

We first fit the space-time ETAS model to the 2005 M5.1 Obsidian Buttes 
earthquake using the Ogata andZhuang [2006] algorithm, in an attempt to account for as 
much of the aftershock behavior as possible. From the Ogata and Zhuang [2006] 
algorithm, theMLE for the ETAS parameters are.fi: = 0.53 events/day/deg^ a=0.92,p = 

1.3, c = 0.01 days, d = 4.8e-5 deg^ q = 2.63, and q = 0.23. However, the data vector 
formed using these parameters resulted in a number of negative seismicity rate values, 
primarily due to the estimate of p, which can lead to instabilities in the filter due to the 
assumption of a Gaussian error distribution. Therefore, we instead use the ordinary-time 
ETAS parameters (.fi: = 0.61 events/day/deg^ a=9M,p = 1.1, and c = 0.001 days). 

The ETAS-predicted rate is subtracted from the observed seismicity rate to form 
the data vector, and the data covariance Rk is estimated. Again because of the tradeoff 
between the parameters rand Aa we fix Aato 1 M Pa, which is roughly consistent with 
faults under hydrostatic conditions at a depth of ~4 km, the depth at which the Obsidian 
Buttes swarm occurred [Chester and Higgs, 1992; Blanpied et al, 1998; Lohman and 
McGuire, 2007]. We then grid-search over values of rand a Using the constraint that 
the integral of the filter estimate of stressing rate must be consistent with the filter 
estimate of y we choose a value of r= 2.5e-4, and maximize the likelihood over values 
of £(Fig. 8). The maximum likelihood estimate occurs when e= 0.5 (Fig. 8a) but again, 

as in the synthetic test, at this value the data covariance multiplier is significantly less 
than 1 (Fig. 8b), and the filter over-fits the data. Our preferred value is £= 0.04, for 

which ~1, and the largest peaks in the data are better fit. 
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4.3 Results 

The stressing rate estinnates for each box using £= 0.04 are shown in Figure 9. 
Figure 10 piots the forward stressing rate estimate for each box on top of one another, 
iiiustrating the fiiter's ability to detect when and in which box the largest transient in the 
region occurs. The largest anomaly occurs in Box 2 and can be associated with a 
geodetically-observed 2005 shallow aseismic creep event on the Obsidian Buttes fault 
that triggered an earthquake swarm [Lohman and McGuire, 2007]. The peak forward 
estimate of stressing rate is 0.042 ±0.004 M Pa/day and the backsmoothed estimate is 
0.022 ±0.006 M Pa/day, roughly two orders of magnitude above tectonic loading. 

The second largest signal also occurred in Box 2 and is related to the Bombay 
Beach earthquake swarm that occurred in M arch 2009. The swarm consisted of ~100s of 
events, the largest of which was a M 4.8 that occurred three days after the swarm initiated. 
Other small anomalies in Boxes 2 and 4 in May 2003 may be related to an earthquake 
swarm that occurred in the Imperial fault zone (located near the boundary of the two 
boxes) [Roland and McGuire , 2009]. Smaller peaks in Box 2 also can be associated with 
earthquake swarms that occurred in the B raw ley seismic zone in 1996, 1998, and 2008 
[Southern California Earthquake Center, http://www.data.scec.org/monthly/index.php]. 

5. Discussion 

Our results highlight the need for a time-dependent background seismicity rate to 
account for variations in seismicity rate due to aseismic processes. Fig. 11 compares the 
observed cumulative number of events in the Salton Trough catalog with the predicted 
number of events from optimizing the space-time ETAS model to the part of the catalog 
that occurred prior to the 2005 Obsidian Buttes swarm and the predicted number of 
events from the filter estimate of seismicity rate. We transform the occurrence times r, of 

the events in the catalog with the theoretical cumulative function t" =|'/l(5)j5, where A 

is the predicted seismicity rate from either ETAS or the filter [Ogata, 1988, 2005]. A 
plot of the cumulative number of events vs. transformed time should be linear if the 
seismicity in the catalog is well described by a particular model. Positive (or negative) 
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deviations from this trend indicate that the modei is under-predicting (or over-predicting) 
the amount of seismicity. Fig. 11 shows that with the ETAS modei, a significant positive 
deviation from this trend occurs near the beginning of the Obsidian Buttes swarm, 
suggesting that anomaious seismicity is occurring that cannot be expiained by the ETAS 
modei aione. Thefiiter prediction however is abie to match the observed seismicity weii. 
Therefore, the seismicity rate anomaiies that appear with respect to the space-time ETAS 
modei, which utiiizes a time-independent background seismicity rate, can be accounted 
for by the time-dependent background seismicity rate produced by our fiiter aigorithm. 

To validate the estimates of S obtained from these seismicity rate variations, we 
can compare our peak stressing rate estimate for the 2005 Obsidian Buttes aseismic 
transient to an estimate based on a geodetically-derived slip model of the deformation 
inverted from InSAR data [Lohman and McGuire, 2007]. The seismicity triggered by 
this transient occurred primarily in the depth range of 4-6 km. We calculate the Coulomb 
stress change [King et al, 1994; Lin and Stein, 2004; Toda et al, 2005] at this depth 
range due to aseismic slip on the shallow part of the fault (Fig. 12). The average total 
Coulomb stress change on this part of the fault is 0.6 M Pa. Based on GPS line-length 
change data, the transient lasted ~1-10 days [Lohman and McGuire, 2007]. Using this 
range of durations, the average stressing rate during the transient then becomes -0.06-0.6 
M Pa/day. For a duration of 5 days (which appears to best describe the GPS data), the 
average stressing rate is 0.12 M Pa/day. From our filter estimate, we obtain a peak 
stressing rate of -0.04 M Pa/day (in the forward filter) and -0.02 M Pa/day (in the 
backsmoothed filter) (Fig. 9). Our estimates from inverting the seismicity catalog are 
within a factor of 5 of the average stressing rate derived from the geodetic data. If we 
take the duration of the transient to match the time step in our filter (20 days) then the 
estimates agree extremely well (0.03 M Pa/day vs. 0.02-0.04 M Pa/day). Given that 
stressing rate increases are likely to be many orders of magnitude over background plate 
tectonic rates [Thatcher, 2001; Toda et al., 2002; Lohman and McGuire, 2007], the 
Salton Trough example demonstrates the feasibility of utilizing our approach to both 
detect and constrain the magnitude of stressing rate transients. 
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6. Conclusion 

We have developed a technique to detect aseismic transients in time and space 
from earthquake catalog data by combining the ETAS and rate-state models of seismicity 
rate into a single data assimilation algorithm to invert catalogs for stressing rate 
variations. We applied it to a synthetic test and successfully detected an aseismic 
transient that involved a factor of 10 increase in stressing rate. We then applied it to a 
catalog from the Salton Trough in California, and successfully detected the onset and 
constrained the absolute magnitude of the largest aseismic transient in a 20 year catalog 
to within a factor of five of the stressing rate estimated with geodetic data. 

Overall, the synthetic test and Salton Trough results suggest that our algorithm is 
a feasible way to detect aseismic stressing rate transients strictly from seismicity catalog 
data. This method may ultimately enable aseismic transient detection in regions lacking 
good geodetic data resolution, such as the (offshore) updip part of subduction zone faults, 
and in time periods prior to the widespread availability of geodetic data. Additionally, a 
seismicity based approach may be more sensitive to small (M4-5) and/or shallow slow- 
slip transients that are not detected by even dense geodetic networks such as the Plate 
Boundary Observatory [Wei et al, 2009]. The results suggest that our seismicity 
inversion method provides an accurate way to detect and locate transient deformation 
strictly from seismicity catalogs and can constrain the absolute magnitude of stressing 
rate variations. 
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Time (days) 


Figure 1. a) Synthetic catalog locations and boxes used in analysis. Parameters used to 
generate the catalog are.fi: = 4.5e-4 events/day/deg^ < 2 = 0.8,jo = 1.2, c = 0.0001 days, d = 
0.001 deg^ <? = 3, r] = 0.3, r = 10 events/day, and Aa= 0.001 M Pa. b) Synthetic catalog 
earthquake magnitudes and occurrence times. Anomalous seismicity triggered by the 
transient is not readily apparent in either the map or the magnitude-time history, which 
are dominated by the aftershock sequences of the largest events in the catalog. 
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Figure 2. a) True stressing rate S used in each spatiai box for the synthetic test. A 
transient rate is imposed in Box 5 (purpie), invoiving an increase in 5 by a factor of 10 
over the background rate of 0.1 M Pa/yr imposed in the other boxes (green), b) Observed 
seismicity rates in each spatiai box (coiors) obtained by binning occurrence times in time 
bins of 1 day. The anomaious seismicity rate due to the transient in Box 5 (purpie) is 
overshadowed by jumps in rate due to aftershock sequences, therefore the stressing rate 
estimate from a straightforward rate-state inversion wiii primariiy refiect changes due to 
these iarger peaks, c) Data vectors (i.e., the aseismicaiiy-triggered seismicity ratei?A) in 
each box obtained by subtracting out the ETAS-predicted rate. The anomaious seismicity 
in Box 5 (purpie) due to the transient is now much more visibie. 
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Figure 3. a) Log likelihood vs. smoothing parameter a b) Data covariance multiplier 
vs. a Ideally should be close to 1, because the data covariance Rk is estimated a 
priori. The maximum likelihood estimate of a\s indicated by a square; our preferred 
value for a, for which ~1, is indicated by a star. 
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Figure 4. Comparison of backsmoothed filter results in Box 5 for different values of s. 
a) Estimates for the rate-state seismicity state variable yfor the M LE, £= 0.79 (teal), and 
our preferred value of £■= 0.16 (black). Also shown are the estimates for )/directly 
derived from the filter (solid lines) and estimates for ^integrated from the filter estimate 
of S using the rate-state model (Eg. 5-6) (dashed lines). Ideally the two estimates at a 
given value of f'should agree with each other, which is better achieved when £= 0.16. b) 
Fit to the observed data vector in Box 5 (red) by each value of £ For the M LE (teal), the 
filter clearly over-fits the data. The £= 0.16 case (black) better fits the largest peak, c) 
Forward filter estimates for S in Box 5 with £= 0.79 (teal) and £= 0.16 (black), 
compared to the true s (red). The £= 0.16 estimate comes closer to recovering the true 
stressing rate. 
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Figure 5. Filter estimates of stressing rate S in each box. T rue stressing rate is shown by 
the red line, the forward estimate by the purple line, and the backsmoothed estimate by 
the black line. Both the forward and backsmoothed results clearly identify in which box 
the transient is located, and produce estimates of the peak stressing rate of the transient 
that are within a factor of 2-3 of the actual peak. 
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Figure 6. Overlay of backsmoothed filter estimates of S in each box (solid colors). The 
true stressing rates are shown by the dashed lines. This emphasizes that the filter is able 
to correctly identify the box in which the transient occurs and estimates the peak stressing 
rate to within a factor of 2-3 of the true peak. 
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Figure?, a) M ap of the Salton Trough region in California showing M > 1.5 seismicity 
occurring from February 1990-August 2009, obtained from the Southern California 
Earthquake Data Center. For our analysis, the region is divided up into the 4 boxes 
indicated, b) M agnitude-time history of the Salton Trough catalog. 
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Figure 8. a) Loglikelihood vs. f'for the Salton Trough data set. b) Data covariance 
multiplier vs. £ Again the maximum likelihood estimate of £(square) differs from 
our preferred value for ^(star), for which As in the synthetic test, the M LE value 
of £(0.5) overfits the data, while the preferred value of £(0.04) fits the largest peaks 
better and also results in better agreement between the direct and integrated filter 
esti mates of y 
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Figure 9. Filter estimates of S in each spatial box for the Salton Trough. The forward 
estimate is shown by the purple line, and the backsmoothed estimate by the black line. 
The largest signal, in Box 2, can be associated with a geodetically-observed aseismic 
transient in the Obsidian Buttes in 2005 [Lohman and McGuire, 2007]. The next largest 
signal, also in Box 2, can be associated with the 2009 Bombay Beach earthquake swarm. 
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Figure 10. Overlay of forward filter estimates of stressing rate in each box (solid colors). 
This emphasizes that the filter is able to detect the largest geodetically-observed transient 
that has occurred in the region (the Obsidian Buttes aseismic creep event in 2005). Other 
smaller anomalies may be related to an earthquake swarm on the Imperial Fault in 2003 
[Roland and McGuire, 2009] in Boxes 2 and 4, earthquake swarms in the B raw ley 
seismic zone in 1996, 1998, and 2008 in Box 2, and the 2009 Bombay Beach swarm in 
Box 2. 
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Figure 11. Cumulative number of events vs. transformed time (i.e., predicted cumulative 
number of events). The red line is a one-to-one line indicating a perfect fit to the 
observed data. ETAS transformed times are calculated using seismicity rates estimated 
from the space-time ETAS model optimized to just prior to the 2005 Obsidian Buttes 
earthquake swarm (event 3779) and extrapolated for the remainder of the catalog (blue 
line). Transformed times are also calculated using the filter estimate of seismicity rate 
(black line). The significant deviation of the blue line from the data shows that the ETAS 
model (with its time-independent background rate) underpredicts the amount of 
seismicity, particularly during the 2005 swarm. The filter estimate (with a time- 
dependent background rate) provides a better fit to the observed cumulative number of 
events. 
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Figure 12. a) Slip model of shallow aseismic creep on the Obsidian Buttes fault inverted 
from InSAR data [Lohman and McGuire, 2007], Black box indicates depth at which 
microseismicity triggered by the transient occurred, b) Calculated Coulomb stress change 
on the fault due to the shallow aseismic slip. The Coulomb stress change averaged over 
the outlined boxes is used to calculate an average stressing rate for the transient (-0.06- 
0.6 M Pa/day for transient durations of 1-10 days). 
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Chapter 5: Detecting A seismic T ransients in the Hokkaido Corner 


Abstract 

Transient aseismic processes alter the stress state of a region and can cause 
seismicity rate anomalies in space and time relative to commonly used aftershock 
models. The presence of such anomalies in subduction zones can therefore indicate stress 
changes occurring due to processes such as fluid flow, afterslip or slow earthquakes. 
Therefore, mapping these anomalies can lead to a better understanding of where and how 
stress accumulates on the megathrust. We have developed a method to invert seismicity 
catalogs for stressing rate variations in time and space. We apply this technique to a 
catalog of events occurring in the Hokkaido corner from 1985-2009. This area consists 
of asperities that rupture in great earthquakes such as the 2003 M8.0 Tokachi-oki 
earthquake. Our initial results identify seismicity anomalies that can be related to 
geodetically-observed afterslip following major interplate earthquakes in the catalog, as 
well as a possible precursory transient before the 1994 M 7.6 earthquake. 
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1. Introduction 

Transient aseismic deformation in subduction zones, such as afterslip following 
major interplate earthquakes or slow slip events on the plate interface, alter the stress 
state of a region and can trigger seismicity [e.g., Hsu et al., 2006; Ozawa ef a/., 2007; 
Delahaye et al., 2009]. This seismicity often does not decay in time or space as quickly 
as typical aftershock sequences, and so appear as seismicity rate anomalies relative to 
models, such as the Epidemic Type Aftershock Sequence (ETAS) model [Ogata, 1988], 
that are based on the empirical laws describing aftershock behavior (e.g., Omori's law, 
Bath's law, and the Gutenberg-Richter magnitude frequency distribution) [L/enos et al., 
2009]. The presence of these seismicity rate anomalies in subduction zones can therefore 
indicate that stress changes are occurring due to processes such as afterslip or slow slip 
events, and identifying regions in which they occur can lead to a better understanding of 
stress accumulation on the megathrust, which has important implications for seismic 
hazard assessment. M oreover, detecting these anomalies could provide a way to map the 
frictional properties of the plate interface by identifying velocity-strengthening regions 
where stress is being released stably through aseismic fault slip as opposed to velocity¬ 
weakening regions where stress is released through recurrent major earthquakes. 

The Hokkaido corner in northeastern japan, where the Pacific plate subducts 
along the Kurile and japan trenches at rates of ~8-9 cm/yr [DeMefs etal., 1990], consists 
of several fixed asperities that rupture in great earthquakes such as the 1994 M7.6 
Sanriku-oki earthquake and the 2003 M 8.0 Tokachi-oki earthquake [e.g., Yamanaka and 
KIkuchl, 2004]. The areas of the plate interface surrounding the asperities release stress 
primarily through aseismic slip [Yamanaka and KIkuchl, 2004; Miyazaki et al., 2004a]. 
A large quantity of both seismic and geodetic data is available in japan, making it a good 
region to investigate seismicity rate anomalies and their relationship to aseismic 
deformation. 

In Chapter 4, we developed a method that can detect seismicity rate anomalies in 
an earthquake catalog and map them to underlying stressing rate variations. We 
implemented a seismicity rate model [Llenos et al., 2009] that combined both the ETAS 
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model and the rate- and state-dependent friction model [Dieterich, 1994] into a data 
assimilation algorithm that inverts seismicity catalogs for stressing rate variations, 
resulting in space-time estimates of the state variables in our model (background stressing 
rate, aseismic stressing rate, and rate-state variable We now apply this technique to a 
catalog of events occurring in the Hokkaido corner from 1985-2009 to detect seismicity 
anomalies and explore how they are related to spatial and temporal variations in stress on 
the megathrust. 


2. M ethod 

We apply the method described in Chapter 4 to analyze the catalog of events from 
northeastern Japan for aseismic stressing rate variations. We assume that the seismicity 
rate R observed in a catalog can be approximated as a linear combination of a component 
Ra (reflecting seismicity triggered by aseismic processes, including long-term tectonic 
loading) and a coseismically-triggered component fic (representing aftershock sequences) 
The component Rc is estimated with the stochastic space-time ETAS model 
[Ogata, 1998; Ogata andZhuang, 2006]. The ETAS model is a point process model 
that represents the seismicity rate R{t,x,y) as a summation of the Omori-like aftershock 
sequences produced by each event prior to time t plus a time-independent background 
seismicity rate jL/[Ogata, 1998; Ogata andZhuang, 2006]: 

R [t, X, y) =iLj{x, y) + ^ /(-(M ,.)- (f -f,.)■ f (x -x,., y -y,; M ,.) (1) 

/ <f 
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The summation in Eq. 1 representing the aftershock sequences in the catalog is 
essentially Rc- The function /iM) is the expected number of events triggered by an 
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earthquake of magnitude M given an aftershock productivity K and scaling parameter a, 
the function ^t) is the probability density function form of the modified Omori law and 
is specified by the Omori decay parameters c and p [Omori, 1894; Utsu, 1961], and the 
function f{x,y;M) describes the spatial distribution of the events triggered by an 
earthquake of magnitude M, specified by the parameters d, q, and scaling parameter q. 

We divide the catalog into a number of spatial boxes, bin the occurrence times in 
each box to obtain seismicity rates in time and space, and subtract the ETAS-predicted Rc 
from the observed rate R in each box to obtain a first-order estimate of a time-dependent 
background seismicity rate Ra that contains variations due to aseismic processes. This 
seismicity rate can be related to a stressing rate on a population of faults in the region 
through the use of the rate- and state-dependent friction model [Dieterich, 1994]: 


r^r 

dr=-^i).-s) (6) 

Aa 

where r is a steady-state reference seismicity rate associated with a reference stressing 
rate 5,, S is the modified Coulomb stressing rate, yls a state variable, and A is a fault 
constitutive parameter. We assume the normal stress a and therefore the frictional 
parameter A aremai n constant. 

We incorporate the rate-state equations into an underlying state-space model. The 
time-dependent state variables used to describe the system are background stressing rate 
Sp, aseismic stressing rate S^,, and the rate-state variable y'm each spatial box. A 

number of time-independent model parameters also must be optimized. Seven of these 
parameters come from the space-time ETAS model, and are optimized using the 
algorithm of Ogata and Zhuang [2006]. Two of the parameters (r and Adj come from 
the rate-state model equations and are estimated a priori. The final two parameters (rand 
£) are the scaling parameters for the random walk processes used to model aseismic 
stressing rate and j/i'espectively and are optimized using a Gaussian likelihood function 
[Harvey, 1989; Segall and Matthews, 1997] subject to the constraints discussed in 
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Chapter 4. We solve the system using an extended Kalman filter algorithm [e.g., Gelb, 
1974; Segall and Matthews, 1997; McGuire and Segall, 2003], which results in time 
histories of the stressing rate in each spatial box. 

3. Data analysis and results 

We apply our method to a catalog of M >3earthquakes that occurred in 
northeastern Japan from 1985-2009, obtained from the japan Meteorological Agency 
(Fig. 1). During this time period, four major interplate earthquakes occurred, three of 
which were off Sanriku (1989 M 7.1, 1992 M 6.9, and 1994 M 7.6) and one of which was 
off Tokachi (2003 M8.0). We divide the region up into 6 spatial boxes and bin the 
occurrence times into time windows of 30 days to obtain the seismicity rate in each box. 
The length of the time step is necessary in order for enough seismicity to occur to allow 
the background rate to be resolved in most of the time steps. We set the rate-state 
frictional parameter A a to be 0.4 M Pa, in agreement with the estimate of Fukuda et al. 
[2009] inferred from GPS data following the Tokachi-oki earthquake, and, following 
Toda and Matsumara [2006] and Ghimire et al. [2009], assume a value for the reference 
stressing rate S, of 0.2 M Pa/yr, which corresponds to a stress drop of 10 M Pa with a 50 
year recurrence interval under full plate coupling [Yamanaka and Kikuchi, 2004; 
Miyazaki etai, 2004a; Yagi, 2004]. 

3.1 ETAS fitting to aftershock sequences 

To explore how well the ETAS model parameters describe the seismicity in the 
catalog, we first fit the ordinary and space-time ETAS models [Ogata, 1988; Ogata and 
Zhuang, 2006] to the aftershock sequences of the four major interplate earthquakes. The 
results suggest that either variations in background rate or in the ETAS parameters over 
time, space and from sequence to sequence could be significant and cause the sequences 
to appear as anomalies from the rest of the catalog. 

Table 1 summarizes the space-time ETAS model parameters that best fit the first 
100 days of aftershocks of the four major interplate earthquakes. These results indicate 
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that the ETAS parameters can vary both in space and in time. The parameters for the first 
two events (1989 M 7.1 and 1992 M 6.9), which occur in the southern part of the catalog, 
are relatively similar, but the parameters for the 1994 M7.6 and 2003 M8.0 are quite 
different. While the background seismicity rate appears to be relatively constant in 
time and space, this may be due to the general insensitivity of the ETAS model to 
background seismicity when fitting individual aftershock sequences [Ogata, 2006]. In 
fact, yt/can vary by several orders of magnitude in both space and time [Ogata et al., 
2003; HainzI and Ogata, 2005; Llenos etal., 2009]. 

For both the 1994 and 2003 events, the ETAS model cannot fit both the early (~5- 
10 days) aftershocks and the late aftershocks (200 days) with the same set of parameters 
(Table 2). We further demonstrate this for theTokachi-oki earthquake by optimizing the 
ordinary ETAS model to the first 5 days of aftershocks and extrapolating the fit to a year¬ 
long catalog of aftershocks. We transform the occurrence times f,- of the events using the 

theoretical cumulative function t" =|'/l(s)ds, where/I is the seismicity rate predicted by 

the optimized ETAS parameters [Ogata, 1988, 2005]. A plot of the cumulative number 
of events in the year-long catalog versus transformed time (i.e., theoretical cumulative 
number) should be linear if the ETAS parameters optimized from the first 5 days fit the 
catalog well. Fig. 2 demonstrates that this is not the case for the 2003 Tokachi-oki 
earthquake. Significantly less earthquakes occur in the later part of the aftershock 
sequence than is predicted by the fit to the earlier part, which suggests that the temporal 
decay in seismicity following the mainshock cannot be well explained by a single set of 
ETAS parameters. 

The space-time ETAS model also has difficulty modeling the spatial decay in 
seismicity in the first month following the 2003 Tokachi-oki event. The ETAS model 
utilizes a spatial probability density function that predicts a power-law decay in 
mainshock-aftershock distance (Eq. 4). However, a histogram of the distances of the 
Tokachi-oki aftershocks reveals that the ETA S model that best fits this sequence does not 
match the observed decay of seismicity in space (Fig. 3). The number of aftershocks in 
fact peaks at distances of ~0.2-0.3 deg in the first month after the mainshock. This 


100 



roughly corresponds with the location of the peak of the 30-day afterslip of this event, 
modeled from GPS data [M/yaza/c; eta/., 2004a] (Fig. 1). 

Taken together, the results of fitting both the ordinary and space-time ETAS 
models to various space-time windows of the catalog show that a single set of ETAS 
parameters alone cannot explain the space-time distribution of seismicity following the 
largest earthquakes in the catalog. Seismicity rate anomalies may arise due to the 
inability of a given set of parameters to model either the spatial or temporal decay in 
aftershocks. This suggests that other triggering mechanisms besides earthquake 
interactions may be responsible for these variations. 

3.2 Filter results 

To form the data vector for the filter, we fit the space-time ETAS model to the 
100-day aftershock sequence of the 2003 Tokachi-oki earthquake in order to account for 
as much of the aftershocks following this event as possible. We then subtract this ETAS- 
predicted rate from the observed rate to estimate the aseismically-triggered seismicity 
rate Ra- This residual thus represents a time-dependent background seismicity rate, and 
becomes the input of the extended Kalman filter. The filter results in estimates of 
background stressing rate Sp, aseismic stressing rate 5 ^, and rate-state variable k As a 

check, the estimate of 5^ can be forward propagated using the rate-state equations (Eq. 
5-6) to obtain the filter's prediction of seismicity rate Ra- This quantity can be combined 
with the ETAS-predicted Rc and integrated to obtain a filter prediction of the cumulative 
number of events in the catalog over time, which can be compared to the observed 
cumulative number of events (Fig. 4). While the fiIter under-predicts the total number of 
events (possibly due to smoothing), it accounts for more of the catalog than the ETAS 
model optimized to the entire catalog alone. This demonstrates that the inclusion of a 
time-dependent background seismicity rate {Ra) is necessary to account for the 
cumulative number of events in the catalog. 

The filter estimates of total stressing rate S (the sum of and Sp) for each box 
are shown in Fig. 5a. The largest signals are clearly due to afterslip following the four 
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interplate earthquakes, which has been observed geodetically. However, two smaller 
anomalies are also detected that were not observed geodetically. In the following 
sections, we will explore these signals in greater detail. 

3.2.1 Afterslip of the 1989 M 7.1 and 1992 M6.9 Sanriku-oki earthquakes 

The 2 November 1989 M7.1 and the 18 July 1992 M6.9 thrust earthquakes 
occurred off the coast of Sanriku in northern Honshu (Fig. 1). The first two large signals 
detected by the filter appear to be spatially and temporally coincident with the afterslip 
associated with these events (Fig. 5b). Although these events occurred prior to the 
establishment of theGEONET continuous GPS network in japan, significant postseismic 
deformation following each was observed in extensometer data [Miura et al., 1993; 
Kawasaki et al., 1995, 2001]. Despite the 1989 earthquake having a larger magnitude 
than the 1992 earthquake, slightly more aseismic moment release occurred for the second 
event [Kawasaki et al., 1995, 2001], which agrees well with the relative sizes of the 
stressing rate peaks we obtain for these two events. 

3.2.2 Afterslip of the 1994 M 7.6 Sanriku-oki earthquake 

A M7.6 earthquake occurred off Sanriku on 12 December 1994 (Fig. 1), with 
rupture propagating westward and the majority of the seismic moment release occurring 
downdip of the epicenter [e.g., Sato et al., 1996; Yagi et al., 2003; Llenos and McGuire, 
2007]. Significant afterslip following this earthquake was detected from both GPS and 
extensometer data [Heki et al., 1997; Kawasaki et al., 2001]. Postseismic slip models 
suggest that the afterslip occurred both in the rupture zone and downdip [Nishimura et 
al., 2000; Yagi et al., 2003]. 

The filter detects large stressing rate anomalies in both Box 4 (downdip of the 
epicenter) and Box 5 (Fig. 5c) that occur following the 1994 earthquake. These transients 
involve peak stressing rates that are ~l-2 orders of magnitude above the tectonic rate. 
Contrary to many of the afterslip models of this event, the transient in Box 5 has a larger 
amplitude and longer duration than the transient detected in Box 4. 
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A particularly interesting anomaly is detected in Box 6, located updip of the 
epicenter near the trench. An increase in stressing rate begins near the start of 1994 that 
peaks ~l-2 months prior to the occurrence of the 1994 earthquake (Fig. 5c). This peak 
stressing rate is a factor of ~10 above the tectonic rate. Compared to the other anomalies 
detected by the filter, this signal is quite small, but represents a large jump in stressing 
rate for the box it is observed in. While not observed geodetically, perhaps because of 
low resolution so far offshore, pre-seismic slip in this region beginning ~8 months prior 
to the 1994 Sanriku-oki earthquake has been inferred by repeating earthquake analysis 
[Uchida et al., 2004]. Repeating earthquakes are assumed to rupture the same small 
asperity, which is surrounded by stable sliding areas; therefore the cumulative slip 
released by the repeating events reflect the aseismic slip release in the surrounding 
regions [A/adeat/ and Johnson, 1998; Igarashi et al., 2003]. The onset of the increased 
stressing rate detected by our filter agrees remarkably well with the accelerated slip 
detected through repeating earthquake analysis. 

3.2.3 Afterslip of the 2003 M8.0 Tokachi-oki earthquake 

A M8.0 earthquake occurred offshore of Hokkaido on 25 September 2003 (Fig. 
1). The abundance of high quality seismic and geodetic data for that event have led to the 
development of detailed coseismic and postseismic slip models [e.g., Yamanaka and 
Kikuchi, 2003; Miyazaki et al., 2004a,b; Yagi, 2004]. The majority of the seismic 
moment release occurred downdip to the northwest of the earthquake's hypocenter. 
Afterslip occurred primarily on the parts of the fault surrounding the coseismic rupture, 
both updip and along-strike [Miyazaki et al., 2004a; Baba et al., 2006]. While the 
afterslip updip of the rupture is not particularly well resolved in land-based geodetic 
models, its occurrence is supported by offshore pressure gauge data [Baba et al., 2006] 
and repeating earthquake analysis [Mafst/bara etal., 2005]. 

The filter detects stressing rate transients in Boxes 2 and 3 (Fig. 5d). However, 
the transient is much larger in Box 3, where the majority of the afterslip occurred, than in 
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Box 2 where the majority of the coseismic slip occurred (Fig. 1). This agrees well with 
thegeodetically-derived afterslip models of this event. 

3.2.4 Afterslip following moderate earthquakes 

Besides the large anomalies associated with afterslip of the major interplate thrust 
events, we have also detected two smaller anomalies in Box 5 that correspond with two 
smaller thrust events: a M 6.5 earthquake that occurred on 8 April 1994 and a M 6.3 that 
occurred on 31 M ay 1998 (Fig. 5b). Afterslip following these events has not been found 
by geodetic data analyses, however repeating earthquake analyses detected accelerations 
in quasi-static slip following both events [Uchida etal., 2003]. The 8 April 1994 event is 
particularly interesting, in that it was located between the hypocenter of the M7.6 that 
occurred later that year (downdip to the west) and the region near the trench (updip to the 
east) where we detected a possible precursory rate change that began in -january- 
February and repeating earthquake analysis detected pre-seismic slow slip [Uchida etal., 
2004]. This supports the idea that aseismic transient deformation initiated near the trench 
in the beginning of 1994 and migrated downdip over the course of the year, culminating 
in a stress concentration around the asperity of the M 7.6 that led to the rupture of that 
asperity in December of that year [Uchida etal., 2004]. 

4. Discussion 

Our results suggest that transient aseismic deformation in subduction zones can be 
detected solely from earthquake catalog data. Besides the examples of afterslip we have 
detected in northeastern japan, recurrent Mw~6.5 slow slip events offshore of the Boso 
peninsula in central japan have also triggered seismicity anomalies [Ozawa et al., 2007; 
Llenos et al., 2009]. Other examples of deformation in subduction zones driving 
seismicity are afterslip following the 2005 Mw 8.7 N las earthquake [Hsu etal., 2006], and 
a slow slip event in the Hikurangi trough in New Zealand that triggered microseismicity 
but no discernible tremor [Delahaye et al., 2009]. These findings clearly suggest that 
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aseismic deformation can trigger seismicity, which may therefore provide a means to 
detect it. 

Our technique may also be useful for mapping spatial variations in frictional 
conditions on the plate interface. For the 2003 Tokachi-oki event, coseismic slip and 
postseismic slip occurred on adjacent patches (Fig. 1), which Miyazaki et ai. [2004a] 
inferred was an indication of spatial heterogeneity in the frictional stability of the plate 
interface. The rupture of the mainshock, which occurs on a velocity-weakening patch, 
loads the adjacent velocity-strengthening patches, which releases the stress stably as 
afterslip. This hypothesis is supported by modeling of GPS data that suggests the timing 
of the early afterslip (~5 hours after the mainshock) is controlled by the stress change due 
to the mainshock and the frictional properties of the megathrust [Fukuda etai., 2009]. 

In our analysis for the Tokachi-oki event, we divided the region such that the 
majority of the coseismic slip occurred in Box 2, while the peak afterslip occurred in Box 
3 (Fig. 1). Although a small transient occurs in Box 2 at the time of the earthquake, the 
transient signal in Box 3 is larger by a factor of ~5 (Fig. 5d). Thus we are clearly able to 
identify the spatial region where the peak afterslip occurred, and distinguish between 
areas that are more velocity-weakening (Box 2) and those that are more velocity¬ 
strengthening (Box 3). With increased spatial resolution (i.e., lower magnitude of 
completeness in the catalog), our technique could provide a way to map the frictional 
conditions on the plate interface simply from earthquake catalog data. 

Finally, we have demonstrated that our method, which relies solely on seismicity 
data, can successfully detect and locate aseismic transients in subduction zones. The 
repeating earthquake analyses off the coast of Hokkaido and Honshu, which also utilize 
seismicity data but are based on entirely different physical models and assumptions, have 
also detected accelerations in quasi-static slip due to aseismic transients [e.g., igarashi et 
ai., 2003; Uchida et ai., 2003, 2004]. These results together suggest that seismicity rate 
histories can provide an alternative way to detect and monitor aseismic deformation in 
areas where land-based geodetic resolution is poor and seafloor geodesy impractical if 
not impossible, such as the offshore updip part of subduction megathrusts. 
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5. Conclusion 

We applied the method developed in Chapter 4 to invert earthquake catalogs for 
stressing rate variations to a catalog of events from northeastern Japan. We successfully 
detected afterslip following four major interplate earthquakes as well as two moderate 
interplate earthquakes. For the 2003 Tokachi-oki event, we were able to spatially 
distinguish between parts of the fault that were velocity-weakening and parts of the fault 
that were velocity-strengthening. We also detected possible precursory aseismic slip 
prior to the 1994 M 7.6 Sanriku-oki earthquake that initiated near the trench and migrated 
downdip, ultimately potentially triggering the rupture of that asperity. Our results 
suggest the following: 1) transient deformation in subduction zones triggers seismicity 
anomalies relative to the ETAS model; 2) our method can successfully detect these 
anomalies in space and time simply from seismic data, which can enable the observation 
of near-trench aseismic deformation; and 3) this method can also be used to map the 
spatial distribution of frictional heterogeneities on the plate interface. 
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Figure 1. Map of the Hokkaido corner showing M >3 seismicity occurring from 1985- 
2009, obtained from the Japan M eteoroiogicai Agency (biack dots). H ypocenters of the 4 
iargest interpiate events are indicated by red stars. Centroids from the Giobai Centroid 
Moment Tensor cataiog for the 1994 M7.6 Sanriku-oki and 2003 M8.0 Tokachi-oki 
earthquakes are indicated by red triangies. A iso shown for these two events are eiiipses 
representing the part of the fauit where the majority of the coseismic moment reiease 
occurred [L/enos and McGuire, 2007]. For the 2003 Tokachi-oki event, a coseismic si ip 
modei (red contours) [Miyazaki etai., 2004b] and 30-day aftersiip modei (biue contours) 
[Miyazaki et a/., 2004a] estimated from GPS are aiso shown. For the fiiter inversion the 
region is divided up into the 6 boxes indicated by the coiored numbers. 
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Figure 2. Results of optimizing the ETAS model to fit the firsts daysof theTokachi-oki 
aftershocks and extrapolated for a year following the earthquake. Bottom panels show the 
magnitudes of the events, a) Observed cumulative number of events in the year following 
the earthquake (blue) compared with the ETAS prediction (red) based on fitting the first 5 
days (dashed line), b) Cumulative number of events versus ETAS transformed time 
(predicted cumulative number of events). The ETAS prediction is shown by the red line 
and observed data in blue. Black lines signify the 2(7bounds of the extrapolation. There 
are significantly less earthquakes in the later part of the sequence than is predicted by the 
parameters that fit the earlier part, suggesting that the same set of parameters cannot be 
used to describe both the early and late parts of the aftershock sequence. 
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Figures. H istogram of mainshock-aftershock distances for aftershocks in Boxes 2 and 3 
in the first month following the 2003 M 8.0 Tokachi-oki earthquake, scaled such that the 
peak occurs at a value of 1. Also shown is the decay predicted by the spatial PDF used in 
the ETAS model (Eg. 4) with parameters estimated for the 2003 earthquake (solid black 
line) and for the 1989 M7.1 event (dashed black line) as a comparison (see Table 1 for 
parameter values). The aftershocks do not decay as quickly with distance as predicted by 
the ETA S model fit to either of the events. 
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Figure 4. Cumulative number of events vs. time observed in the catalog (heavy line), 
compared with the number of events predicted from the ETAS model alone (dashed line), 
and from the filter estimate of stressing rate (thin line). While still under-predicting the 
total number of events, the filter estimate which utilizes a time-dependent background 
seismicity rate does a better job at fitting the observations than the ETAS model which 
utilizes a time-independent background rate. 
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Figure 5. a) Forward filter estimates of stressing rate S in each of the spatial boxes 
indicated in Fig. 1. The largest signals detected are associated with geodetically- 
observed afterslip following major interplate earthquakes in 1989 (M7.1), 1992 (M6.9), 
1994 (M7.6) and 2003 (M8.0). b) Results for Box 5, showing geodetically-observed 
afterslip following the 1989, 1992, and 1994 off-Sanriku earthquakes, as well as afterslip 
following a 1994 M6.5 and 1998 M6.3 event, which has been inferred from repeating 
earthquake analyses [Uchida et al., 2003; Uchida et a/., 2004], c) Filter estimates of 5 
for the updip (Box 6), downdip (Box 4), and hypocentral (Box 5) areas of the 1994 M 7.6 
Sanriku-oki earthquake. Afterslip is concentrated in Box 4-5. A possible preseismic 
transient begins in early 1994 updip of the rupture (Box 6) and peaks ~2 months before 
the rupture, d) Results for the 2003 M8.0 Tokachi-oki event. A larger signal occurs in 
Box 3 (where the peak afterslip occurred [Miyazaki etai., 2004a]) than in Box 2 (where 
thecoseismic slip occurred [Miyazaki etai., 2004b]). 
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